pointgrp.h

//
// 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
 */
 
#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 {
 
// //////////////////////////////////////////////////////////////////
 
class SymmetryOperation {
  private:
    double d[3][3];
 
  public:
    SymmetryOperation();
    SymmetryOperation(const SymmetryOperation &);
    ~SymmetryOperation();
 
    double trace() const { return d[0][0]+d[1][1]+d[2][2]; }
 
    double* operator[](int i) { return d[i]; }
 
    const double* operator[](int i) const { return d[i]; }
 
    double& operator()(int i, int j) { return d[i][j]; }
 
    double operator()(int i, int j) const { return d[i][j]; }
 
    void zero() { memset(d,0,sizeof(double)*9); }
 
    SymmetryOperation operate(const SymmetryOperation& r) const;
 
    SymmetryOperation transform(const SymmetryOperation& r) const;
    
    void unit() { zero(); d[0][0] = d[1][1] = d[2][2] = 1.0; }
 
    void E() { unit(); }
    
    void i() { zero(); d[0][0] = d[1][1] = d[2][2] = -1.0; }
 
    void sigma_h() { unit(); d[2][2] = -1.0; }
 
    void sigma_xz() { unit(); d[1][1] = -1.0; }
 
    void sigma_yz() { unit(); d[0][0] = -1.0; }
 
    void rotation(int n);
    void rotation(double theta);
    
    void c2_x() { i(); d[0][0] = 1.0; }
 
    void c2_y() { i(); d[1][1] = 1.0; }
 
    void transpose();
 
    void print(std::ostream& =ExEnv::out0()) const;
};
 
// //////////////////////////////////////////////////////////////////
 
class SymRep {
  private:
    int n;
    double d[5][5];
 
  public:
    SymRep(int =0);
    SymRep(const SymmetryOperation&);
    ~SymRep();
 
    operator SymmetryOperation() const;
    
    inline double trace() const;
 
    void set_dim(int i) { n=i; }
    
    double* operator[](int i) { return d[i]; }
    const double* operator[](int i) const { return d[i]; }
 
    double& operator()(int i, int j) { return d[i][j]; }
    double operator()(int i, int j) const { return d[i][j]; }
 
    void zero() { memset(d,0,sizeof(double)*25); }
 
    SymRep operate(const SymRep& r) const;
 
    SymRep transform(const SymRep& r) const;
    
    void unit() {
      zero(); d[0][0] = d[1][1] = d[2][2] = d[3][3] = d[4][4] = 1.0;
    }
    
    void E() { unit(); }
    
    void i() { zero(); d[0][0] = d[1][1] = d[2][2] = d[3][3] = d[4][4] = -1.0;}
 
    void sigma_h();
 
    void sigma_xz();
 
    void sigma_yz();
 
    void rotation(int n);
    void rotation(double theta);
    
    void c2_x();
 
    void c2_y();
 
    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;
 
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&);
    IrreducibleRepresentation(int,int,const char*,const char* =0);
 
    ~IrreducibleRepresentation();
 
    IrreducibleRepresentation& operator=(const IrreducibleRepresentation&);
 
    void init(int =0, int =0, const char* =0, const char* =0);
    
    int order() const { return g; }
 
    int degeneracy() const { return degen; }
 
    int complex() const { return complex_; }
 
    int nproj() const { return degen*degen; }
 
    int nrot() const { return nrot_; }
 
    int ntrans() const { return ntrans_; }
 
    const char * symbol() const { return symb; }
 
    const char * symbol_ns() const { return (csymb?csymb:symb); }
 
    double character(int i) const {
      return complex_ ? 0.5*rep[i].trace() : rep[i].trace();
    }
 
    double p(int x1, int x2, int i) const { return rep[i](x1,x2); }
    
    double p(int d, int i) const {
      int dc=d/degen; int dr=d%degen;
      return rep[i](dr,dc);
    }
 
    void print(std::ostream& =ExEnv::out0()) const;
};
 
// ///////////////////////////////////////////////////////////
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
 
    int parse_symbol();
    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();
    CharacterTable(const char*);
    CharacterTable(const char*,const SymmetryOperation&);
 
    CharacterTable(const CharacterTable&);
    ~CharacterTable();
 
    CharacterTable& operator=(const CharacterTable&);
 
    int nirrep() const { return nirrep_; }
    int order() const { return g; }
    const char * symbol() const { return symb; }
    IrreducibleRepresentation& gamma(int i) { return gamma_[i]; }
    SymmetryOperation& symm_operation(int i) { return symop[i]; }
 
    int complex() const {
      if (pg==CN || pg==SN || pg==CNH || pg==T || pg==TH)
        return 1;
      return 0;
    }
 
    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;
    }
 
    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;
    }
 
    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;
    }
    
    void print(std::ostream& =ExEnv::out0()) const;
};
 
// ///////////////////////////////////////////////////////////
 
class PointGroup: public SavableState {
  private:
    std::string symb;
    SymmetryOperation frame;
    SCVector3 origin_;
 
  public:
    PointGroup();
    PointGroup(std::string);
    PointGroup(std::string,SymmetryOperation&);
    PointGroup(std::string,SymmetryOperation&,const SCVector3&);
    PointGroup(const Ref<KeyVal>&);
 
    PointGroup(StateIn&);
    PointGroup(const PointGroup&);
    PointGroup(const Ref<PointGroup>&);
    ~PointGroup();
 
    PointGroup& operator=(const PointGroup&);
 
    int order() const { return char_table().order(); }
 
    bool equiv(const Ref<PointGroup> &, double tol = 1.0e-6) const;
 
    CharacterTable char_table() const;
    std::string symbol() const { return symb; }
    SymmetryOperation& symm_frame() { return frame; }
    const SymmetryOperation& symm_frame() const { return frame; }
    SCVector3& origin() { return origin_; }
    const SCVector3& origin() const { return origin_; }
 
    void set_symbol(std::string);
 
    void save_data_state(StateOut& so);
 
    void print(std::ostream&o=ExEnv::out0()) const;
};
 
}
 
#endif
 
// Local Variables:
// mode: c++
// c-file-style: "ETS"
// End: