maketab.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

//
// maketab.cc
//
// Modifications are
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Edward Seidl <[e-mail]>
// 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: [e-mail]
 *      June, 1993
 */

#include <math.h>
#include <stdio.h>
#include <string.h>

#include <math/symmetry/pointgrp.h>
#include <util/misc/formio.h>

using namespace std;
using namespace sc;

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif


/*
 * 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();