maketab.cc

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

        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: