#include <vector>
#include "string.h"
#include "phonon.h"
#include "green.h"
#include "timer.h"
#include "global.h"

#ifdef UseSPG
extern "C"{
#include "spglib.h"
}
#endif

/*------------------------------------------------------------------------------
 * Private method to evaluate the phonon dispersion curves
 *----------------------------------------------------------------------------*/
void Phonon::pdisp()
{
  // ask the output file name and write the header.
  char str[MAXLINE];
  for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
#ifdef UseSPG
  // ask method to generate q-lines
  int method = 2;
  printf("Please select your method to generate the phonon dispersion:\n");
  printf("  1. Manual, should always work;\n");
  printf("  2. Automatic, works only for 3D crystals (CMS49-299).\nYour choice [2]: ");
  if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f"));
  method = 2 - method%2;
  printf("Your  selection: %d\n", method);
#endif
  printf("\nPlease input the filename to output the dispersion data [pdisp.dat]:");
  if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdisp.dat");
  char *ptr = strtok(str," \t\n\r\f");
  char *fname = new char[strlen(ptr)+1];
  strcpy(fname,ptr);

  // to store the nodes of the dispersion curve
  std::vector<double> nodes;        nodes.clear();
  std::vector<std::string> ndstr;   ndstr.clear();
  std::vector<double *> qs, qe;     qs.clear(); qe.clear();
  std::vector<int> nqbin;           nqbin.clear();

  // now the calculate the dispersion curve
  double qstr[3], qend[3];
  int nq = MAX(MAX(dynmat->nx,dynmat->ny),dynmat->nz)/2+1;
  qend[0] = qend[1] = qend[2] = 0.;

  double *egvs = new double [ndim];
#ifdef UseSPG
  if (method == 1){
#endif
    while (1){
      for (int i = 0; i < 3; ++i) qstr[i] = qend[i];
  
      int quit = 0;
      printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]);
      int n = count_words(fgets(str, MAXLINE, stdin));
      ptr = strtok(str, " \t\n\r\f");
      if ((n == 1) && (strcmp(ptr,"q") == 0)) break;
      else if (n >= 3){
        qstr[0] = atof(ptr);
        qstr[1] = atof(strtok(NULL, " \t\n\r\f"));
        qstr[2] = atof(strtok(NULL, " \t\n\r\f"));
      }
  
      do printf("Please input the end q-point in unit of B1->B3: ");
      while (count_words(fgets(str, MAXLINE, stdin)) < 3);
      qend[0] = atof(strtok(str,  " \t\n\r\f"));
      qend[1] = atof(strtok(NULL, " \t\n\r\f"));
      qend[2] = atof(strtok(NULL, " \t\n\r\f"));
  
      printf("Please input the # of points along the line [%d]: ", nq);
      if (count_words(fgets(str, MAXLINE, stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f"));
      nq = MAX(nq,2);
  
      double *qtmp = new double [3];
      for (int i = 0; i < 3; ++i) qtmp[i] = qstr[i];
      qs.push_back(qtmp);

      qtmp = new double [3];
      for (int i = 0; i < 3; ++i) qtmp[i] = qend[i];
      qe.push_back(qtmp);

      nqbin.push_back(nq);
      ndstr.push_back("");
    }
    ndstr.push_back("");

#ifdef UseSPG
  } else {

    memory->grow(atpos, dynmat->nucell, 3, "pdisp:atpos");
    memory->grow(attyp, dynmat->nucell,    "pdisp:attyp");
  
    num_atom = dynmat->nucell;
    // set default, in case system dimension under study is not 3.
    for (int i = 0; i < dynmat->nucell; ++i)
    for (int idim = 0; idim < 3; ++idim) atpos[i][idim] = 0.;
    for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j) latvec[i][j] = 0.;
    for (int i = 0; i < 3; ++i) latvec[i][i] = 1.;

    // get atomic type info
    for (int i = 0; i < num_atom; ++i) attyp[i] = dynmat->attyp[i];

    // get unit cell vector info
    int ndim = 0;
    for (int idim = 0; idim < 3; ++idim)
    for (int jdim = 0; jdim < 3; ++jdim) latvec[jdim][idim] = dynmat->basevec[ndim++];

    // get atom position in unit cell; fractional
    for (int i = 0; i < num_atom; ++i)
    for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim];

    // display the unit cell info read
    printf("\n");
    for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n");
    printf("The basis vectors of the unit cell:\n");
    for (int idim = 0; idim < 3; ++idim) printf("  A%d = %lg %lg %lg\n", idim+1, latvec[0][idim], latvec[1][idim], latvec[2][idim]);
    printf("Atom(s) in the unit cell:\n");
    printf("  No.  type  sx  sy sz\n");
    for (int i = 0; i < MIN(num_atom, NUMATOM); ++i) printf("  %d %d %lg %lg %lg\n", i+1, attyp[i], atpos[i][0], atpos[i][1], atpos[i][2]);
    if (num_atom > NUMATOM) printf("  ... (%d atoms omitted.)\n", num_atom-NUMATOM);

    char symbol[11];
    double symprec = 1.e-4, pos[num_atom][3];
    for (int i = 0; i < num_atom; ++i)
    for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j];

    int spgnum  = spg_get_international(symbol, latvec, pos, attyp, num_atom, symprec);
    printf("The space group number of your unit cell is: %d => %s\n", spgnum, symbol);
    for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n");
  
    // angles
    const double la = sqrt(latvec[0][0] * latvec[0][0] + latvec[0][1] * latvec[0][1] + latvec[0][2] * latvec[0][2]);
    const double lb = sqrt(latvec[1][0] * latvec[1][0] + latvec[1][1] * latvec[1][1] + latvec[1][2] * latvec[1][2]);
    const double lc = sqrt(latvec[2][0] * latvec[2][0] + latvec[2][1] * latvec[2][1] + latvec[2][2] * latvec[2][2]);
    const double cosa = sqrt(latvec[1][0] * latvec[2][0] + latvec[1][1] * latvec[2][1] + latvec[1][2] * latvec[2][2])/(lb*lc);
    const double cosg = sqrt(latvec[0][0] * latvec[1][0] + latvec[0][1] * latvec[1][1] + latvec[0][2] * latvec[1][2])/(la*lb);
    double ivec[3][3];
    ndim = 0;
    for (int idim = 0; idim < 3; ++idim)
    for (int jdim = 0; jdim < 3; ++jdim) ivec[jdim][idim] = dynmat->ibasevec[ndim++];
    const double ka = sqrt(ivec[0][0] * ivec[0][0] + ivec[0][1] * ivec[0][1] + ivec[0][2] * ivec[0][2]);
    const double kb = sqrt(ivec[1][0] * ivec[1][0] + ivec[1][1] * ivec[1][1] + ivec[1][2] * ivec[1][2]);
    const double coskg = sqrt(ivec[0][0] * ivec[1][0] + ivec[0][1] * ivec[1][1] + ivec[0][2] * ivec[1][2])/(ka*kb);
    
    double *qtmp;
    if (spgnum >= 1 && spgnum <= 2){  // Triclinic system
      if (fabs(coskg) > ZERO){        // A.14, TRI1a and TRI2a

        ndstr.push_back("X");
        // X-G
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] =  qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");

        // G-Y
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("Y/L");

        // L-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");

        // G-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("Z/N");

        // N-G
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");

        // G-M
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("M/R");

        // R-G
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = qtmp[1] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");

      } else {        // A.14, TRI1b and TRI2b

        ndstr.push_back("X");
        // X-G
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = -0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] =  qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");

        // G-Y
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[1] = qtmp[2] = 0.; qtmp[0] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("Y/L");

        // L-G
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");

        // G-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = -0.5; qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("Z/N");

        // N-G
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = -0.5; qtmp[1] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");

        // G-M
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("M/R");

        // R-G
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = -0.5; qtmp[2] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");
      }

    } else if (spgnum >= 3 && spgnum <= 15){ // Monoclinic

      if (symbol[0] == 'P'){             // MCL-P
        const double eta = (1.-lb*cosa/lc)/(2.*(1.-cosa*cosa));
        const double niu = 0.5 - eta * lc * cosa / lb;

        ndstr.push_back("{/Symbol G}");
        // G-Y
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("Y");

        // Y-H
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = eta; qtmp[2] = 1.-niu;
        qe.push_back(qtmp);

        ndstr.push_back("H");

        // H-C
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = eta; qtmp[2] = 1.-niu;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("C");

        // C-E
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("E");

        // E-M1
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = 1.-eta; qtmp[2] = niu;
        qe.push_back(qtmp);

        ndstr.push_back("M_1");

        // M1-A
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = 1.-eta; qtmp[2] = niu;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("A");

        // A-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("X");

        // X-H1
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[2] = 1.-eta; qtmp[1] = niu;
        qe.push_back(qtmp);

        ndstr.push_back("H_1/M");

        // M-D
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[2] = eta; qtmp[1] = 1.-niu;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("D");

        // D-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("Z/Y");

        // Y-D
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("D");

      } else {                  // MCL-C

        if (coskg < 0.){        // MCLC1
          const double xi = (2. - lb * cosa / lc) / (4.*(1.-cosa*cosa));
          const double eta = 0.5 + 2.*xi*lc/lb*cosa;
          const double psi = 0.75 - la * la /(4.*lb*lb*(1.-cosa*cosa));
          const double phi = psi + (0.75 - psi) * lb / lc * cosa;

          ndstr.push_back("{/Symbol G}");
          // G-Y
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
  
          ndstr.push_back("Y");
  
          // Y-F
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 1.-xi; qtmp[2] = 1.-eta;
          qe.push_back(qtmp);
  
          ndstr.push_back("F");
  
          // F-L
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 1.-xi; qtmp[2] = 1.-eta;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
          qe.push_back(qtmp);
  
          ndstr.push_back("L");
  
          // L-I
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = phi; qtmp[1] = 1.-phi; qtmp[2] = 0.5;
          qe.push_back(qtmp);
  
          ndstr.push_back("I/I_1");
  
          // I1-Z
          qtmp = new double [3];
          qtmp[0] = 1.-phi; qtmp[1] = -qtmp[0]; qtmp[2] = 0.5;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qe.push_back(qtmp);
  
          ndstr.push_back("Z");
  
          // Z-F1
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = xi; qtmp[2] = eta;
          qe.push_back(qtmp);
  
          ndstr.push_back("F_1/Y");
  
          // Y-X1
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = psi; qtmp[1] = 1.-psi; qtmp[2] = 0.;
          qe.push_back(qtmp);
  
          ndstr.push_back("X_1/X");
  
          // X-G
          qtmp = new double [3];
          qtmp[0] = 1.-psi; qtmp[1] = -qtmp[0]; qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
  
          ndstr.push_back("{/Symbol G}");
  
          // G-N
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
  
          ndstr.push_back("N/M");
  
          // M-G
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
  
          ndstr.push_back("{/Symbol G}");
  
        } else if ( fabs(coskg) < ZERO) {     // MCLC2
          const double xi = (2. - lb * cosa / lc) / (4.*(1.-cosa*cosa));
          const double eta = 0.5 + 2.*xi*lc/lb*cosa;
          const double psi = 0.75 - la * la /(4.*lb*lb*(1.-cosa*cosa));
          const double phi = psi + (0.75 - psi) * lb / lc * cosa;

          ndstr.push_back("{/Symbol G}");
          // G-Y
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
  
          ndstr.push_back("Y");
  
          // Y-F
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 1.-xi; qtmp[2] = 1.-eta;
          qe.push_back(qtmp);
  
          ndstr.push_back("F");
  
          // F-L
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 1.-xi; qtmp[2] = 1.-eta;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
          qe.push_back(qtmp);
  
          ndstr.push_back("L");
  
          // L-I
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = phi; qtmp[1] = 1.-phi; qtmp[2] = 0.5;
          qe.push_back(qtmp);
  
          ndstr.push_back("I/I_1");
  
          // I1-Z
          qtmp = new double [3];
          qtmp[0] = 1.-phi; qtmp[1] = -qtmp[0]; qtmp[2] = 0.5;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qe.push_back(qtmp);
  
          ndstr.push_back("Z");
  
          // Z-F1
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = xi; qtmp[2] = eta;
          qe.push_back(qtmp);
  
          ndstr.push_back("F_1/N");
  
          // N-G
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
  
          ndstr.push_back("{/Symbol G}");
  
          // G-M
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
  
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
  
          ndstr.push_back("M");
  
        } else {

          double flag = lb / lc * cosa + lb * lb / (la * la *(1.-cosa*cosa));
          if (fabs(flag) < ZERO){           // MCLC4
            const double miu = 0.25*(1. + lb * lb / (la *la));
            const double del = lb * lc * cosa / (2.*la*la);
            const double xi  = miu - 0.25 + (1. - lb * cosa / lc)/(4.*(1.-cosa*cosa));
            const double eta = 0.5 + 2.*xi*lc/lb*cosa;
            const double phi = 1. + xi - 2.*miu;
            const double psi = eta - 2.*del;
   
            ndstr.push_back("{/Symbol G}");
            // G-Y
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = miu; qtmp[2] = del;
            qe.push_back(qtmp);
    
            ndstr.push_back("Y");

            // Y-F
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = miu; qtmp[2] = del;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 1.-phi; qtmp[2] = 1.-psi;
            qe.push_back(qtmp);
    
            ndstr.push_back("F");

            // F-H
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 1.-phi; qtmp[2] = 1.-psi;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = xi; qtmp[2] = eta;
            qe.push_back(qtmp);
    
            ndstr.push_back("H");

            // H-Z
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = xi; qtmp[2] = eta;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
            qe.push_back(qtmp);
    
            ndstr.push_back("Z");

            // Z-I
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
            qe.push_back(qtmp);
    
            ndstr.push_back("I/H_1");

            // H1-Y1
            qtmp = new double [3];
            qtmp[0] = 1.-xi; qtmp[1] = -xi; qtmp[2] = 1.-eta;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del;
            qe.push_back(qtmp);
    
            ndstr.push_back("Y_1");

            // Y1-X
            qtmp = new double [3];
            qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("X");

            // X-G
            qtmp = new double [3];
            qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("{/Symbol G}");

            // G-N
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("N/M");

            // M-G
            qtmp = new double [3];
            qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("{/Symbol G}");

          } else if (flag < 1.){            // MCLC3
            const double miu = 0.25*(1. + lb * lb / (la *la));
            const double del = lb * lc * cosa / (2.*la*la);
            const double xi  = miu - 0.25 + (1. - lb * cosa / lc)/(4.*(1.-cosa*cosa));
            const double eta = 0.5 + 2.*xi*lc/lb*cosa;
            const double phi = 1. + xi - 2.*miu;
            const double psi = eta - 2.*del;
   
            ndstr.push_back("{/Symbol G}");
            // G-Y
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = miu; qtmp[2] = del;
            qe.push_back(qtmp);
    
            ndstr.push_back("Y");

            // Y-F
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = miu; qtmp[2] = del;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 1.-phi; qtmp[2] = 1.-psi;
            qe.push_back(qtmp);
    
            ndstr.push_back("F");

            // F-H
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 1.-phi; qtmp[2] = 1.-psi;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = xi; qtmp[2] = eta;
            qe.push_back(qtmp);
    
            ndstr.push_back("H");

            // H-Z
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = xi; qtmp[2] = eta;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
            qe.push_back(qtmp);
    
            ndstr.push_back("Z");

            // Z-I
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
            qe.push_back(qtmp);
    
            ndstr.push_back("I/H_1");

            // I-F1
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = phi; qtmp[2] = phi - 1.; qtmp[1] = psi;
            qe.push_back(qtmp);
    
            ndstr.push_back("F_1/H_1");

            // H1-Y1
            qtmp = new double [3];
            qtmp[0] = 1.-xi; qtmp[1] = -xi; qtmp[2] = 1.-eta;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del;
            qe.push_back(qtmp);
    
            ndstr.push_back("Y_1");

            // Y1-X
            qtmp = new double [3];
            qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("X");

            // X-G
            qtmp = new double [3];
            qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("{/Symbol G}");

            // G-N
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("N/M");

            // M-G
            qtmp = new double [3];
            qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("{/Symbol G}");

          } else {                          // MCLC5
            const double xi = (lb*lb/(la*la) + (1.-lb*cosa/lc)/(1.-cosa*cosa))*0.25;
            const double eta = 0.5 + 2.*xi*lc*cosa/lb;
            const double miu = 0.5*eta + lb * lb /(4.*la*la) - lb*lc/(2.*la*la)*cosa;
            const double niu = 2.*miu - xi;
            const double omg = (4.*niu - 1. - lb*lb*(1.-cosa*cosa)/(la*la))*lc/(2.*lb*cosa);
            const double del = xi*lc*cosa/lb + omg*0.5 - 0.25;
            const double rho = 1.-xi*la*la/(lb*lb);
   
            ndstr.push_back("{/Symbol G}");
            // G-Y
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = miu; qtmp[2] = del;
            qe.push_back(qtmp);
    
            ndstr.push_back("Y");

            // Y-F
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = miu; qtmp[2] = del;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = niu; qtmp[2] = omg;
            qe.push_back(qtmp);
    
            ndstr.push_back("F");

            // F-L
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = niu; qtmp[2] = omg;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
            qe.push_back(qtmp);
    
            ndstr.push_back("Y");

            // L-I
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = rho; qtmp[1] = 1.-rho; qtmp[2] = 0.5;
            qe.push_back(qtmp);
    
            ndstr.push_back("I/I_1");

            // I1-Z
            qtmp = new double [3];
            qtmp[0] = 1.-rho; qtmp[1] = -qtmp[0]; qtmp[2] = 0.5;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
            qe.push_back(qtmp);
    
            ndstr.push_back("Z");

            // Z-H
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = xi; qtmp[2] = eta;
            qe.push_back(qtmp);
    
            ndstr.push_back("H");

            // H-F1
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = xi; qtmp[2] = eta;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 1.-niu; qtmp[2] = 1.-omg;
            qe.push_back(qtmp);
    
            ndstr.push_back("F_1/H_1");

            // H1-Y1
            qtmp = new double [3];
            qtmp[0] = 1.-xi; qtmp[1] = -xi; qtmp[2] = 1.-eta;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del;
            qe.push_back(qtmp);
    
            ndstr.push_back("Y_1");

            // Y1-X
            qtmp = new double [3];
            qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("Y_1");

            // X-G
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("{/Symbol G}");

            // G-N
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("N/M");

            // M-G
            qtmp = new double [3];
            qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
            qs.push_back(qtmp);
    
            qtmp = new double [3];
            qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
            qe.push_back(qtmp);
    
            ndstr.push_back("{/Symbol G}");
          }
        }
      }

    } else if (spgnum >= 16 && spgnum <= 74){ // Orthorhombic

      if (symbol[0] == 'P'){             // ORC

        ndstr.push_back("{/Symbol G}");
        // G-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("X");

        // X-S
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("S");

        // S-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("{/Symbol G}");

        // G-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("Z");

        // Z-U
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("U");

        // U-R
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("R");

        // R-T
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("T");

        // T-Z
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("Z/Y");

        // Y-T
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("T/U");

        // U-X
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("X/S");

        // S-R
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("R");

      } else if (symbol[0] == 'F'){        // ORCF

        double flag = 1./(la*la) - 1./(lb*lb) - 1./(lc*lc);
        if (fabs(flag) < ZERO){            // ORCF3

          const double xi  = 0.25 * (1. + la*la/(lb*lb) - la*la/(lc*lc));
          const double eta = 0.25 * (1. + la*la/(lb*lb) + la*la/(lc*lc));

          ndstr.push_back("{/Symbol G}");
          // G-Y
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Y");

          // Y-T
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5;
          qe.push_back(qtmp);
      
          ndstr.push_back("T");

          // T-Z
          qtmp = new double [3];
          qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Z");

          // Z-G
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("{/Symbol G}");

          // G-X
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta;
          qe.push_back(qtmp);
      
          ndstr.push_back("X");

          // X-A1
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5-xi; qtmp[2] = 1.-xi;
          qe.push_back(qtmp);
      
          ndstr.push_back("A_1");

          // A1-Y
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5-xi; qtmp[2] = 1.-xi;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Y/X");

          // X-A
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5+xi; qtmp[2] = xi;
          qe.push_back(qtmp);
      
          ndstr.push_back("A");

          // A-Z
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5+xi; qtmp[2] = xi;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Z/L");

          // L-G
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("{/Symbol G}");

        } else if (flag > 0.){             // ORCF1

          const double xi  = 0.25 * (1. + la*la/(lb*lb) - la*la/(lc*lc));
          const double eta = 0.25 * (1. + la*la/(lb*lb) + la*la/(lc*lc));

          ndstr.push_back("{/Symbol G}");
          // G-Y
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Y");

          // Y-T
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5;
          qe.push_back(qtmp);
      
          ndstr.push_back("T");

          // T-Z
          qtmp = new double [3];
          qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Z");

          // Z-G
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("{/Symbol G}");

          // G-X
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta;
          qe.push_back(qtmp);
      
          ndstr.push_back("X");

          // X-A1
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5-xi; qtmp[2] = 1.-xi;
          qe.push_back(qtmp);
      
          ndstr.push_back("A_1");

          // A1-Y
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5-xi; qtmp[2] = 1.-xi;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Y/T");

          // T-X1
          qtmp = new double [3];
          qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 1.; qtmp[1] = qtmp[2] = 1.-eta;
          qe.push_back(qtmp);
      
          ndstr.push_back("X_1/X");

          // X-A
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5+xi; qtmp[2] = xi;
          qe.push_back(qtmp);
      
          ndstr.push_back("A");

          // A-Z
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5+xi; qtmp[2] = xi;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Z/L");

          // L-G
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("{/Symbol G}");

        } else {                           // ORCF2

          const double eta = 0.25 * (1. + la*la/(lb*lb) - la*la/(lc*lc));
          const double phi = 0.25 * (1. + lc*lc/(lb*lb) - lc*lc/(la*la));
          const double del = 0.25 * (1. + lb*lb/(la*la) - lb*lb/(lc*lc));

          ndstr.push_back("{/Symbol G}");
          // G-Y
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Y");

          // Y-C
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[2] = 0.5-eta; qtmp[1] = 1.-eta;
          qe.push_back(qtmp);
      
          ndstr.push_back("C");

          // C-D
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[2] = 0.5-eta; qtmp[1] = 1.-eta;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.5-del; qtmp[2] = 0.5; qtmp[1] = 1.-del;
          qe.push_back(qtmp);
      
          ndstr.push_back("D");

          // D-X
          qtmp = new double [3];
          qtmp[0] = 0.5-del; qtmp[2] = 0.5; qtmp[1] = 1.-del;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
          qe.push_back(qtmp);
      
          ndstr.push_back("X");

          // X-G
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("{/Symbol G}");

          // G-Z
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Z");

          // Z-D1
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.5+del; qtmp[1] = 0.5; qtmp[2] = del;
          qe.push_back(qtmp);
      
          ndstr.push_back("D_1");

          // D1-H
          qtmp = new double [3];
          qtmp[0] = 0.5+del; qtmp[1] = 0.5; qtmp[2] = del;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 1.-phi; qtmp[1] = 0.5-phi; qtmp[2] = 0.5;
          qe.push_back(qtmp);
      
          ndstr.push_back("H");

          // H-C
          qtmp = new double [3];
          qtmp[0] = 1.-phi; qtmp[1] = 0.5-phi; qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5-eta; qtmp[2] = 1.-eta;
          qe.push_back(qtmp);
      
          ndstr.push_back("C/C_1");

          // C1-Z
          qtmp = new double [3];
          qtmp[0] = 0.5; qtmp[1] = 0.5+eta; qtmp[2] = eta;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Z/X");

          // X-H1
          qtmp = new double [3];
          qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = phi; qtmp[1] = 0.5+phi; qtmp[2] = 0.5;
          qe.push_back(qtmp);
      
          ndstr.push_back("H_1/H");

          // H-Y
          qtmp = new double [3];
          qtmp[0] = 1.-phi; qtmp[1] = 0.5-phi; qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("Y/L");

          // L-G
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
      
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
      
          ndstr.push_back("{/Symbol G}");
        }

      } else if (symbol[0] == 'C'){        // ORCC

        const double xi = 0.25 * (1. + la*la/(lb*lb));

        ndstr.push_back("{/Symbol G}");
        // G-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = xi; qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("X");

        // X-S
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = xi; qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("S");

        // S-R
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("R");

        // R-A
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = xi; qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("A");

        // A-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = xi; qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("Z");

        // Z-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("{/Symbol G}");

        // G-Y
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = -0.5; qtmp[1] = 0.5; qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("Y");

        // Y-X1
        qtmp = new double [3];
        qtmp[0] = -0.5; qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = -xi; qtmp[1] = 1.-xi; qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("X_1");

        // X1-A1
        qtmp = new double [3];
        qtmp[0] = -xi; qtmp[1] = 1.-xi; qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = -xi; qtmp[1] = 1.-xi; qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("A_1");

        // A1-T
        qtmp = new double [3];
        qtmp[0] = -xi; qtmp[1] = 1.-xi; qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("T");

        // T-Y
        qtmp = new double [3];
        qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = -0.5; qtmp[1] = 0.5; qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("Y/Z");

        // Z-T
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("T");

      } else {                             // ORCI

        const double xi  = 0.25 * (1. + la*la/(lc*lc));
        const double eta = 0.25 * (1. + lb*lb/(lc*lc));
        const double del = (lb*lb-la*la)/(4.*lc*lc);
        const double miu = (lb*lb+la*la)/(4.*lc*lc);

        ndstr.push_back("{/Symbol G}");
        // G-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = -xi; qtmp[1] = qtmp[2] = xi;
        qe.push_back(qtmp);
    
        ndstr.push_back("X");

        // X-L
        qtmp = new double [3];
        qtmp[0] = -xi; qtmp[1] = qtmp[2] = xi;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = -miu; qtmp[1] = miu; qtmp[2] = 0.5-del;
        qe.push_back(qtmp);
    
        ndstr.push_back("L");

        // L-T
        qtmp = new double [3];
        qtmp[0] = -miu; qtmp[1] = miu; qtmp[2] = 0.5-del;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("T");

        // T-W
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
        qe.push_back(qtmp);
    
        ndstr.push_back("W");

        // W-R
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("R");

        // R-X1
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = xi; qtmp[1] = 1.-xi; qtmp[2] = -xi;
        qe.push_back(qtmp);
    
        ndstr.push_back("X_1");

        // X1-Z
        qtmp = new double [3];
        qtmp[0] = xi; qtmp[1] = 1.-xi; qtmp[2] = -xi;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("Z");

        // Z-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("{/Symbol G}");

        // G-Y
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = eta; qtmp[1] = -eta;
        qe.push_back(qtmp);
    
        ndstr.push_back("Y");

        // Y-S
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = eta; qtmp[1] = -eta;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
    
        ndstr.push_back("S");

        // S-W
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
        qe.push_back(qtmp);
    
        ndstr.push_back("W/L_1");

        // L1-Y
        qtmp = new double [3];
        qtmp[0] = miu; qtmp[1] = -miu; qtmp[2] = 0.5+del;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = eta; qtmp[1] = -eta;
        qe.push_back(qtmp);
    
        ndstr.push_back("Y/Y_1");

        // Y1-Z
        qtmp = new double [3];
        qtmp[0] = 1.-eta; qtmp[1] = eta; qtmp[2] = -eta;
        qs.push_back(qtmp);
    
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5;
        qe.push_back(qtmp);
    
        ndstr.push_back("Z");
      }

    } else if (spgnum >= 75 && spgnum <= 142){ // Tetragonal

      if (symbol[0] == 'P'){             // TET

        ndstr.push_back("{/Symbol G}");
        // G-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("X");

        // X-M
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("M");

        // M-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);

        ndstr.push_back("{/Symbol G}");

        // G-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("Z");

        // Z-R
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("R");

        // R-A
        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("A");

        // A-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("Z/X");

        // X-R
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("R/M");

        // M-A
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);

        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);

        ndstr.push_back("A");

      } else {                           // BCT
        if (la > lc){                    // BCT1
          const double eta = 0.25 * (1. + lc*lc/(la*la));

          ndstr.push_back("{/Symbol G}");
          // G-X
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("X");

          // X-M
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("M");

          // M-G
          qtmp = new double [3];
          qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
   
          ndstr.push_back("{/Symbol G}");

          // G-Z
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = eta; qtmp[2] = -eta;
          qe.push_back(qtmp);
   
          ndstr.push_back("Z");

          // Z-P
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = eta; qtmp[2] = -eta;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
          qe.push_back(qtmp);
   
          ndstr.push_back("P");

          // P-N
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("N");

          // N-Z1
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = -eta; qtmp[1] = 1.-eta; qtmp[2] = eta;
          qe.push_back(qtmp);
   
          ndstr.push_back("Z_1");

          // Z1-M
          qtmp = new double [3];
          qtmp[0] = -eta; qtmp[1] = 1.-eta; qtmp[2] = eta;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("M");

          // X-P
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
          qe.push_back(qtmp);
   
          ndstr.push_back("P");

        } else {                         // BCT2

          const double eta = 0.25 * (1. + la*la/(lc*lc));
          const double xi  = la*la/(2.*lc*lc);

          ndstr.push_back("{/Symbol G}");
          // G-X
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("X");

          // X-Y
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = -xi; qtmp[1] = xi; qtmp[2] = 0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("Y");

          // Y-Sigma
          qtmp = new double [3];
          qtmp[0] = -xi; qtmp[1] = xi; qtmp[2] = 0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = -eta; qtmp[1] = qtmp[2] = eta;
          qe.push_back(qtmp);
   
          ndstr.push_back("{/Symbol S}");

          // Sigma-G
          qtmp = new double [3];
          qtmp[0] = -eta; qtmp[1] = qtmp[2] = eta;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qe.push_back(qtmp);
   
          ndstr.push_back("{/Symbol G}");

          // G-Z
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = qtmp[2] = 0.;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("Z");

          // Z-Sigma_1
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = eta; qtmp[1] = 1.-eta; qtmp[2] = -eta;
          qe.push_back(qtmp);
   
          ndstr.push_back("{/Symbol S}_1");

          // Sigma_1-N
          qtmp = new double [3];
          qtmp[0] = eta; qtmp[1] = 1.-eta; qtmp[2] = -eta;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("N");

          // N-P
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = qtmp[1] = 0.25;
          qe.push_back(qtmp);
   
          ndstr.push_back("P");

          // P-Y1
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = qtmp[1] = 0.25;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -xi;
          qe.push_back(qtmp);
   
          ndstr.push_back("Y_1");

          // Y1-Z
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -xi;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
          qe.push_back(qtmp);
   
          ndstr.push_back("Z");

          // X-P
          qtmp = new double [3];
          qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
          qs.push_back(qtmp);
   
          qtmp = new double [3];
          qtmp[0] = qtmp[2] = qtmp[1] = 0.25;
          qe.push_back(qtmp);
   
          ndstr.push_back("P");
        }
      }
    } else if (spgnum >= 143 && spgnum <= 167){ // Trigonal
      if (cosg > 0.){                           // RHL1

        const double eta = (1.+4.*cosa)/(2.+4.*cosa);
        const double niu = 0.75 - 0.5*eta;

        ndstr.push_back("{/Symbol G}");
        // G-L
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
 
        ndstr.push_back("L");

        // L-B1
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = 1.-eta; qtmp[2] = eta - 1.;
        qe.push_back(qtmp);
 
        ndstr.push_back("B_1/B");

        // B-Z
        qtmp = new double [3];
        qtmp[0] = eta; qtmp[1] = 0.5; qtmp[2] = 1.-eta;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
 
        ndstr.push_back("Z");

        // Z-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
 
        ndstr.push_back("{/Symbol G}");

        // G-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = niu; qtmp[1] = 0.; qtmp[2] = -niu;
        qe.push_back(qtmp);
 
        ndstr.push_back("X/Q");

        // Q-F
        qtmp = new double [3];
        qtmp[0] = 1.-niu; qtmp[1] = niu; qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qe.push_back(qtmp);
 
        ndstr.push_back("F");

        // F-P1
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 1.-niu; qtmp[2] = 1.-eta;
        qe.push_back(qtmp);
 
        ndstr.push_back("P_1");

        // P1-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 1.-niu; qtmp[2] = 1.-eta;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
 
        ndstr.push_back("Z/L");

        // L-P
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = eta; qtmp[1] = qtmp[2] = niu;
        qe.push_back(qtmp);
 
        ndstr.push_back("P");

      } else {                                  // RHL2
        const double eta = 0.5*(1.+cosa)/(1.-cosa);
        const double niu = 0.75 - 0.5*eta;

        ndstr.push_back("{/Symbol G}");
        // G-P
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 1.-niu; qtmp[1] = -niu;
        qe.push_back(qtmp);
 
        ndstr.push_back("P");

        // P-Z
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 1.-niu; qtmp[1] = -niu;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
        qe.push_back(qtmp);
 
        ndstr.push_back("Z");

        // Z-Q
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = eta;
        qe.push_back(qtmp);
 
        ndstr.push_back("Q");

        // Q-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = eta;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
 
        ndstr.push_back("{/Symbol G}");

        // G-F
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.;
        qe.push_back(qtmp);
 
        ndstr.push_back("F");

        // F-P1
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = niu; qtmp[1] = qtmp[2] = niu - 1.;
        qe.push_back(qtmp);
 
        ndstr.push_back("P_1");

        // P1-Q1
        qtmp = new double [3];
        qtmp[0] = niu; qtmp[1] = qtmp[2] = niu - 1.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = 1.-eta; qtmp[1] = qtmp[2] = -eta;
        qe.push_back(qtmp);
 
        ndstr.push_back("Q_1");

        // Q1-L
        qtmp = new double [3];
        qtmp[0] = 1.-eta; qtmp[1] = qtmp[2] = -eta;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
 
        ndstr.push_back("L");

        // L-Z
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
 
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
        qe.push_back(qtmp);
 
        ndstr.push_back("Z");
      }

    } else if (spgnum >= 168 && spgnum <= 194){ // Hexagonal

      ndstr.push_back("{/Symbol G}");
      // G-M
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = qtmp[2] = 0.;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
      qe.push_back(qtmp);
 
      ndstr.push_back("M");
     
      // M-K
      qtmp = new double [3];
      qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.;
      qe.push_back(qtmp);
 
      ndstr.push_back("K");
     
      // K-G
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = qtmp[2] = 0.;
      qe.push_back(qtmp);
 
      ndstr.push_back("{/Symbol G}");
     
      // G-A
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = qtmp[2] = 0.;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
      qe.push_back(qtmp);
 
      ndstr.push_back("A");
     
      // A-L
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
      qe.push_back(qtmp);
 
      ndstr.push_back("L");
     
      // L-H
      qtmp = new double [3];
      qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.5;
      qe.push_back(qtmp);
 
      ndstr.push_back("H");
     
      // H-A
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.5;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
      qe.push_back(qtmp);
 
      ndstr.push_back("A/L");
     
      // L-M
      qtmp = new double [3];
      qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.;
      qe.push_back(qtmp);
 
      ndstr.push_back("M/K");
     
      // K-H
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.;
      qs.push_back(qtmp);
 
      qtmp = new double [3];
      qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.5;
      qe.push_back(qtmp);
 
      ndstr.push_back("H");
     
    } else if (spgnum >= 195 && spgnum <= 230){ // Cubic

      if (symbol[0] == 'P'){                    // CUB

        ndstr.push_back("{/Symbol G}");
        // G-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("X");
     
        // X-M
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qe.push_back(qtmp);
   
        ndstr.push_back("M");
     
        // M-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
   
        ndstr.push_back("{/Symbol G}");
     
        // G-R
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("R");
     
        // R-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("X/M");
     
        // M-R
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("R");

      } else if (symbol[0] == 'F'){          // FCC

        ndstr.push_back("{/Symbol G}");
        // G-X
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qe.push_back(qtmp);
   
        ndstr.push_back("X");
     
        // X-W
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = 0.25; qtmp[2] = 0.75;
        qe.push_back(qtmp);
   
        ndstr.push_back("W");
     
        // W-K
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = 0.25; qtmp[2] = 0.75;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.375; qtmp[2] = 0.75;
        qe.push_back(qtmp);
   
        ndstr.push_back("K");
     
        // K-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.375; qtmp[2] = 0.75;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
   
        ndstr.push_back("{/Symbol G}");
     
        // G-L
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("L");
     
        // L-U
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = 0.625; qtmp[1] = 0.25; qtmp[2] = 0.625;
        qe.push_back(qtmp);
   
        ndstr.push_back("U");
     
        // U-W
        qtmp = new double [3];
        qtmp[0] = 0.625; qtmp[1] = 0.25; qtmp[2] = 0.625;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = 0.25; qtmp[2] = 0.75;
        qe.push_back(qtmp);
   
        ndstr.push_back("W");
     
        // W-L
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = 0.25; qtmp[2] = 0.75;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("L");
     
        // L-K
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.5;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.375; qtmp[2] = 0.75;
        qe.push_back(qtmp);
   
        ndstr.push_back("K/U");
     
        // U-X
        qtmp = new double [3];
        qtmp[0] = 0.625; qtmp[1] = 0.25; qtmp[2] = 0.625;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = 0.5; qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("X");
     
      } else {                               // BCC

        ndstr.push_back("{/Symbol G}");
        // G-H
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("H");
     
        // H-N
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("N");
     
        // N-G
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qe.push_back(qtmp);
   
        ndstr.push_back("{/Symbol G}");
     
        // G-P
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
        qe.push_back(qtmp);
   
        ndstr.push_back("P");
     
        // P-H
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("H/P");
     
        // P-N
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = qtmp[2] = 0.25;
        qs.push_back(qtmp);
   
        qtmp = new double [3];
        qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5;
        qe.push_back(qtmp);
   
        ndstr.push_back("N");
     
      }

    } else {
      printf("\nSorry, failed to identify the crystal system, please use the manual mode.\n");
    }

    // to determine the number of points along each line, with a step size of 0.05
    const double qs_inv = 1./QSTEP;
    int nbin = qs.size();
    if (nbin > 0) printf("\nPhonon dispersion will be evaluated along lines:\n\t%s", ndstr[0].c_str());
    for (int is = 0; is < nbin; ++is){
      double *qstr = qs[is];
      double *qend = qe[is];
      double ql = 0.;
      for (int i = 0; i < 3; ++i) ql += (qend[i] - qstr[i])*(qend[i] - qstr[i]);
      int nqpt = MAX(int(sqrt(ql) * qs_inv + 0.5), 2);
      nqbin.push_back(nqpt);

      printf("-%s", ndstr[is+1].c_str());
    }
    if (nbin > 0) printf("\n");
  }
#endif

  FILE *fp = fopen(fname, "w");
  fprintf(fp,"# q     qr    freq\n");
  fprintf(fp,"# 2pi/L  2pi/L %s\n", dynmat->funit);

  double qr = 0., dq, q[3], qinc[3];
  int nbin = qs.size();
  for (int is = 0; is < nbin; ++is){
    double *qstr = qs[is];
    double *qend = qe[is];
    int nbin = nqbin[is];
    for (int i = 0; i < 3; ++i) qinc[i] = (qend[i]-qstr[i])/double(nbin-1);
    dq = sqrt(qinc[0]*qinc[0]+qinc[1]*qinc[1]+qinc[2]*qinc[2]);
  
    nodes.push_back(qr);
    for (int i = 0; i < 3; ++i) q[i] = qstr[i];
    for (int ii = 0; ii < nbin; ++ii){
      double wii = 1.;
      dynmat->getDMq(q, &wii);
      if (wii > 0.){
        dynmat->geteigen(egvs, 0);
        fprintf(fp,"%lg %lg %lg %lg ", q[0], q[1], q[2], qr);
        for (int i = 0; i < ndim; ++i) fprintf(fp," %lg", egvs[i]);
      }
      fprintf(fp,"\n");
  
      for (int i = 0; i < 3; ++i) q[i] += qinc[i];
      qr += dq;
    }
    qr -= dq;
    delete []qstr;
    delete []qend;
  }
  qs.clear(); qe.clear();
  if (qr > 0.) nodes.push_back(qr);
  fclose(fp);
  delete []egvs;

  // write the gnuplot script which helps to visualize the result
  int nnd = nodes.size();
  if (nnd > 1){
    const char qmk = char(34); // "
    fp = fopen("pdisp.gnuplot", "w");
    fprintf(fp,"set term post enha colo 20\nset out %cpdisp.eps%c\n\n", qmk, qmk);
    fprintf(fp,"set xlabel %cq%c\n", qmk, qmk);
    fprintf(fp,"set ylabel %cfrequency (THz)%c\n\n", qmk, qmk);
    fprintf(fp,"set xrange [0:%lg]\nset yrange [0:*]\n\n", nodes[nnd-1]);
    fprintf(fp,"set grid xtics\n");
    fprintf(fp,"# {/Symbol G} will give you letter gamma in the label\nset xtics (");
    for (int i = 0; i < nnd-1; ++i) fprintf(fp,"%c%s%c %lg, ", qmk, ndstr[i].c_str(), qmk, nodes[i]);
    fprintf(fp, "%c%s%c %lg)\n\n", qmk, ndstr[nnd-1].c_str(), qmk, nodes[nnd-1]);
    fprintf(fp, "unset key\n\n");
    fprintf(fp, "plot %c%s%c u 4:5 w l lt 1", qmk, fname, qmk);
    for (int i = 1; i < ndim; ++i) fprintf(fp,",\\\n%c%c u 4:%d w l lt 1", qmk, qmk, i+5);
    fclose(fp);

    printf("\nPhonon dispersion data are written to: %s, you can visualize the results\n", fname);
    printf("by invoking: `gnuplot pdisp.gnuplot; gv pdisp.eps`\n");
  }
  for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");

  delete []fname;
  nodes.clear();
  ndstr.clear();

return;
}