/**********************************************************************
  analysis_example.c:

     analysis_example.c provides examples for analyzing and 
     unitilizing Kohn-Sham Hamiltonian, overlap, and density
     matrices which are stored in filename.scfout.

  Log of read_scfout.c:

     2/July/2003  Released by T.Ozaki 

  ******************************************************************
    You can utilize a filename.scfout which is generated by the SCF
    calculation of ABRED by the following procedure:

     1. Define your main routine as follows:

           int main(int argc, char *argv[]) 

     2. Include a header file, "read_scfout.h", in your main
        (if you want, also in other routines) as follows:

           #include "read_scfout.h"
  
     3. Call a function, read_scfout(), in the main routine as follows:

           read_scfout(argv);
  ******************************************************************

***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "read_scfout.h"
#include "Inputtools.h"

typedef long int INTEGER;

void Eigen_lapack(double **a, double *ko, int n);
void Draw_Gcube(double coef[6]);
void read_coordinates(char *argv[]);

int *Choo_atom ;  // array to store the glbal number of choosen atoms

int main(int argc, char *argv[]) 
{

  static int ct_AN,h_AN,Gh_AN,i,j,TNO1,TNO2;  
  static int spin,Rn,num0,num1,num2,num3;

  int II,JJ ;   // dummy variables for loop
  int Anum;     // temporary variable
  int Number_Choo ;  // the number of Choosen atom
  int TNumOrbs,TNumOrbs3;
  int SPI,SPJ;
  int optEV ;  // variable to be used for option of eigenvector-printing
  int gcube_on;
  int *MP ;  // array which specify a head position of a full matrix
  double ***SmallHks, ***O_SmallHks ;  
  double **SmallOLP, **O_SmallOLP ;  

  /* variable & arrays for PART-2; same with that of Cluster_DFT.c */

  static int l,n,n2,n1,i1,j0,j1,k1,l1;
  static int P_min,num_eigen ;
  static double **ko, *M1;
  static double **B, ***C, **D;
  static double sum,sum1,Tsum;                                                       // v
  static double coef[6];
  static double **LDOS;

  read_scfout(argv);
  read_coordinates(argv);

  printf("Number of choosen atoms?\n");
  scanf("%d",&Number_Choo);

  if (Number_Choo<=0){
    printf("Invalid number\n");
  }

  /************************************************************************* 
     PART-1 :  Constructing the selected-full Hamiltonian & Overlap matrix
  **************************************************************************/

  /*
    allocation of arrays:

    int Choo_atom[Number_Choo]; 
    int MP[Number_Choo]; 
    double **SmallHks ;  
    double **SmallOLP ;  
  */

  Choo_atom = (int*)malloc(sizeof(int)*Number_Choo); 
  MP = (int*)malloc(sizeof(int)*Number_Choo); 

  /* read Choo_atom */

  printf("Specify choosen atoms\n");
  for (i=0; i<Number_Choo; i++){
    scanf("%d",&Choo_atom[i]);
  }

  /*
   make an array MP which specify the starting
   position of atom II in the martix such as
   a full but small Hamiltonian
  */

  Anum = 1;
  for (i=0; i<Number_Choo; i++){
    MP[i] = Anum;
    ct_AN = Choo_atom[i];
    Anum = Anum + Total_NumOrbs[ct_AN];    
  }
  TNumOrbs = Anum - 1;
  TNumOrbs3 = TNumOrbs + 3;  

  for (i=0; i<Number_Choo; i++){
    ct_AN = Choo_atom[i];
    printf("i=%i ct_AN=%i Total_NumOrbs=%i MP=%i\n",
           i,ct_AN,Total_NumOrbs[ct_AN],MP[i]);
  }

  /*
    allocation of arrays:

    double **SmallHks ;  
    double **SmallOLP ;  
  */

  SmallHks = (double***)malloc(sizeof(double**)*(SpinP_switch+1)); 
  for (spin=0; spin<=SpinP_switch; spin++){
    SmallHks[spin] = (double**)malloc(sizeof(double*)*TNumOrbs3); 
    for (i=0; i<TNumOrbs3; i++){
      SmallHks[spin][i] = (double*)malloc(sizeof(double)*TNumOrbs3); 
    }
  }

  SmallOLP = (double**)malloc(sizeof(double*)*TNumOrbs3); 
  for (i=0; i<TNumOrbs3; i++){
    SmallOLP[i] = (double*)malloc(sizeof(double)*TNumOrbs3); 
  }

  O_SmallHks = (double***)malloc(sizeof(double**)*(SpinP_switch+1)); 
  for (spin=0; spin<=SpinP_switch; spin++){
    O_SmallHks[spin] = (double**)malloc(sizeof(double*)*TNumOrbs3); 
    for (i=0; i<TNumOrbs3; i++){
      O_SmallHks[spin][i] = (double*)malloc(sizeof(double)*TNumOrbs3); 
    }
  }

  O_SmallOLP = (double**)malloc(sizeof(double*)*TNumOrbs3); 
  for (i=0; i<TNumOrbs3; i++){
    O_SmallOLP[i] = (double*)malloc(sizeof(double)*TNumOrbs3); 
  }

  // sorting ?

  for (spin=0; spin<=SpinP_switch; spin++){
    //    printf("Kohn-Sham Hamiltonian spin=%i\n",spin);              // v

    for (II=0; II<Number_Choo; II++){
      SPI = MP[II];  
      for (JJ=0; JJ<Number_Choo; JJ++){
        SPJ = MP[JJ];  
	ct_AN = Choo_atom[II] ;
	TNO1 = Total_NumOrbs[ct_AN];
	for (h_AN=0; h_AN<=FNAN[ct_AN]; h_AN++){
	  Gh_AN = natn[ct_AN][h_AN];
	  if (Gh_AN == Choo_atom[JJ]){
	    Rn = ncn[ct_AN][h_AN];
	    TNO2 = Total_NumOrbs[Gh_AN];

	    for (i=0; i<TNO1; i++){
	      for (j=0; j<TNO2; j++){
		SmallHks[spin][i+SPI][j+SPJ] = Hks[spin][ct_AN][h_AN][i][j]; 
		SmallOLP[i+SPI][j+SPJ] = OLP[ct_AN][h_AN][i][j]; 
	      }
	    }

	    /*
	    printf("glbal index=%i  local index=%i (grobal=%i, Rn=%i)\n",
		   ct_AN,h_AN,Gh_AN,Rn);
	    for (i=0; i<TNO1; i++){
	      for (j=0; j<TNO2; j++){
		printf("%10.7f ",Hks[spin][ct_AN][h_AN][i][j]); 
	      }
	      printf("\n");
	    }
	    */

	  }
	}
      }
    }
  }

  /* store original selected Hks and S */

  for (spin=0; spin<=SpinP_switch; spin++){
    printf("spin=%i Full Hamiltonian matrix of selected atoms\n",spin);
    for (i=1; i<=TNumOrbs; i++){
      for (j=1; j<=TNumOrbs; j++){
        printf("%7.3f ",SmallHks[spin][i][j]);
        // store original information
	O_SmallHks[spin][i][j] = SmallHks[spin][i][j]; 
      }
      printf("\n");
    }
  }

  for (i=1; i<=TNumOrbs; i++){
    for (j=1; j<=TNumOrbs; j++){
      printf("%7.3f ",SmallOLP[i][j]);
      O_SmallOLP[i][j] = SmallOLP[i][j];   // store original information
    }
    printf("\n");
  }

  for (spin=0; spin<=SpinP_switch; spin++){
    printf("spin=%i Full Hamiltonian matrix of selected atoms\n",spin);
    for (i=1; i<=TNumOrbs; i++){
      for (j=1; j<=TNumOrbs; j++){
        printf("%7.3f ",SmallHks[spin][i][j]);
      }
      printf("\n");
    }
  }

  printf("Full overlap matrix of selected atoms\n");
  for (i=1; i<=TNumOrbs; i++){
    for (j=1; j<=TNumOrbs; j++){
      printf("%7.3f ",SmallOLP[i][j]);
    }
    printf("\n");
  }

  /************************************************************************* 

     PART-2 :  Starting the part that diagonalize the selected-full 

                                                Hamiltonian & Overlap matrix

  **************************************************************************/

  /*******************************************
   allocation of arrays:

   double ko[SpinP_switch+1][TNumOrbs3];
   double M1[TNumOrbs3];
   double B[TNumOrbs3][TNumOrbs3];

  ********************************************/

  ko = (double**)malloc(sizeof(double*)*(SpinP_switch+1));
  for (spin=0; spin<=SpinP_switch; spin++){
    ko[spin] = (double*)malloc(sizeof(double)*TNumOrbs3); 
  }
  M1 = (double*)malloc(sizeof(double)*TNumOrbs3); 

  B = (double**)malloc(sizeof(double*)*TNumOrbs3); 
  for (i=0; i<TNumOrbs3; i++){
    B[i] = (double*)malloc(sizeof(double)*TNumOrbs3); 
  }

  C = (double***)malloc(sizeof(double**)*(SpinP_switch+1)); 
  for (spin=0; spin<=SpinP_switch; spin++){
    C[spin] = (double**)malloc(sizeof(double*)*TNumOrbs3); 
    for (i=0; i<TNumOrbs3; i++){
      C[spin][i] = (double*)malloc(sizeof(double)*TNumOrbs3); 
    }
  }

  D = (double**)malloc(sizeof(double*)*TNumOrbs3); 
  for (i=0; i<TNumOrbs3; i++){
    D[i] = (double*)malloc(sizeof(double)*TNumOrbs3); 
  }

  /*******************************************
   diagonalize the overlap matrix

   first 
   SmallOLP -> OLP matrix
   after call Eigen_lapack
   SmallOLP -> eigenvectors of OLP matrix
  ********************************************/

  Eigen_lapack(SmallOLP,ko[0],TNumOrbs);
  for (l=1; l<=TNumOrbs; l++){
    M1[l] = 1.0/sqrt(ko[0][l]);
  }

  /*
    for (l=1; l<=TNumOrbs; l++){
    printf("%i %15.12f\n",l,M1[l]);
    }
  */


  /****************************************************
    Calculations of eigenvalues for up and down spins
  ****************************************************/

  n = TNumOrbs;

  for (spin=0; spin<=SpinP_switch; spin++){
    for (i1=1; i1<=n; i1++){
      for (j1=1; j1<=n; j1++){
        sum = 0.0;
        for (l=1; l<=n; l++){
	  sum = sum + SmallHks[spin][i1][l]*SmallOLP[l][j1]*M1[j1]; 
        }
        C[spin][i1][j1] = sum;
      }
    }

    for (i1=1; i1<=n; i1++){
      for (j1=1; j1<=n; j1++){
        sum = 0.0;
        for (l=1; l<=n; l++){
	  sum = sum + M1[i1]*SmallOLP[l][i1]*C[spin][l][j1];
	  //sum = sum + M1[i1]*SmallOLP[l][i1]*C[spin][l][j1];
        }
        B[i1][j1] = sum;
      }
    }

    for (i1=1; i1<=n; i1++){
      for (j1=1; j1<=n; j1++){
        D[i1][j1] = B[i1][j1];    
      }
    }

    Eigen_lapack(D,ko[spin],n);

    /****************************************************
        Transformation to the original eigen vectors.
                          NOTE 244P
    ****************************************************/

    for (i1=1; i1<=n; i1++){
      for (j1=1; j1<=n; j1++){
        C[spin][i1][j1] = 0.0;
      }
    }

    for (i1=1; i1<=n; i1++){
      for (j1=1; j1<=n; j1++){
        sum = 0.0;
        for (l=P_min; l<=n; l++){
          sum = sum + SmallOLP[i1][l]*M1[l]*D[l][j1];       
        }
        C[spin][i1][j1] = sum;
      }
    }
  }

  for (spin=0; spin<=SpinP_switch; spin++){
    printf("\nspin=%i \n",spin);                                       // v  
    for (i=1; i<=TNumOrbs; i++){
      printf("%ith eigenvalue of HC=eSC: %15.12f\n",i,ko[spin][i]);
    }
  }

  /*
  for (spin=0; spin<=SpinP_switch; spin++){
    printf("C spin=%i\n",spin);
    for (i1=1; i1<=n; i1++){
      for (j1=1; j1<=n; j1++){
        printf("%7.4f ",C[spin][i1][j1]);
      }
      printf("\n");
    }
  }
  */

  /****************************************************
                Part for Checking 'HC=eSC'
  ****************************************************/

  for (spin=0; spin<=SpinP_switch; spin++){
    printf("spin=%i \n", spin);
    for (j1=1; j1<=n; j1++){
      for (i1=1; i1<=n; i1++){
	sum = 0.0;
	sum1= 0.0;
	for (l=1; l<=n; l++){
          sum = sum + O_SmallHks[spin][i1][l]*C[spin][l][j1];
	  sum1 = sum1 + O_SmallOLP[i1][l]*C[spin][l][j1];    // *ko[spin][j1];
	}
	sum1 = ko[spin][j1] * sum1 ;    // ko[spin][i1]??
        Tsum = Tsum + fabs(sum-sum1);
      }
      printf("Check ko=%i |HC-eSC|=%15.12f\n",j1,Tsum);
    }
  }

  /*
  for (spin=0; spin<=SpinP_switch; spin++){
    printf("spin=%i \n", spin);
    for (i1=1; i1<=n; i1++){
      for (j1=1; j1<=n; j1++){
	sum = 0.0;
	sum1= 0.0;
	for (l=1; l<=n; l++){
	sum = sum + O_SmallHks[spin][i1][l]*C[spin][l][j1];
	sum1 = sum1 + O_SmallOLP[i1][l]*C[spin][l][j1];    // *ko[spin][j1];
	}
	sum1 = ko[spin][j1] * sum1 ;    // ko[spin][i1]??
	printf("l.h.s - r.h.s = %10.7f\n", sum-sum1 );
      }
    }
  }
  */

  /****************************************************
              printing out the eigenvectors
  ****************************************************/

  printf("\nDo you want eigenvectors also? (yes:1 / no:0)");
  scanf("%i",&optEV);
  if (optEV == 1){
    num0 = 7; 
    for (spin=0; spin<=SpinP_switch; spin++){
      printf("\nspin=%i \n",spin);                                       // v  

      num1 = TNumOrbs/num0;
      num2 = TNumOrbs%num0;

      for (i1=0; i1<=num1; i1++){
        j0 = i1*num0;
        for (i=-2; i<=TNumOrbs; i++){
          for (j=0; j<=num0; j++){
            j1 = j0 + j;
            if (j1<=TNumOrbs){
              if (i==-2){
                if (j==0)
                  printf("    ");
                else
                  printf("%7d     ",j1);
	      }
              else if (i==-1){
                if (j==0)
                  printf("     ");
                else
                  printf("%10.7f  ",ko[spin][j1]);
	      }
              else if (i==0){ 
                if (j==0)
                  printf("    ");
                else 
                 printf("        ");
	      }
              else {
                if (j==0)
                  printf("%4d ",i);
                else 
                  printf("%10.7f  ",C[spin][i][j1]);
	      }
	    }
	  }
          printf("\n");
	}
      } 

    }
  }

  /*
  printf("\nDo you want eigenvectors also? (yes:1 / no:0)");
  scanf("%i",&optEV);
  if (optEV == 1){
    for (spin=0; spin<=SpinP_switch; spin++){
      printf("\nspin=%i \n",spin);                                       // v  
      for (i=1; i<=TNumOrbs; i++){
	printf("%ith eigenvector: ",i);
	printf("{");
	for (j=1; j<=TNumOrbs; j++){
	  printf("%15.12f,",C[spin][i][j]);
	}
	printf("}\n");
      }
    }
  }
  */

  /****************************************************
       making Gaussian cube data of MO(d-orbitals) 
  ****************************************************/

  printf("\nDo you want Gcube of MO (d-orbitals)? (yes:1 / no:0)");
  scanf("%i",&gcube_on);
  if (gcube_on == 1){

    printf("\nWhich eigenstate? (yes:1 / no:0)");
    scanf("%i",&num_eigen);
    printf("\n up or down? (up:0 / down:0)");
    scanf("%i",&spin);

    coef[1] = C[spin][10][num_eigen];
    coef[2] = C[spin][11][num_eigen];
    coef[3] = C[spin][12][num_eigen];
    coef[4] = C[spin][13][num_eigen];
    coef[5] = C[spin][14][num_eigen];
    Draw_Gcube(coef);
  }

  /****************************************************
               calculate PDOS of Mn atom
  ****************************************************/

  LDOS = (double**)malloc(sizeof(double*)*30); 
  for (i=0; i<30; i++){
    LDOS[i] = (double*)malloc(sizeof(double)*TNumOrbs3); 
  }

  for (spin=0; spin<=SpinP_switch; spin++){
    printf("\nPDOS spin=%i\n",spin);
    for (i1=10; i1<=14; i1++){
      for (i=1; i<=TNumOrbs; i++){
        LDOS[i1][i] = 0.0;
        for (j=1; j<=TNumOrbs; j++){
          LDOS[i1][i] += C[spin][i1][i]*C[spin][j][i]*O_SmallOLP[j][i1]; 
	}
      }
    }

    for (i=1; i<=TNumOrbs; i++){
      sum = 0.0; 
      for (i1=10; i1<=14; i1++) sum = sum + LDOS[i1][i];
      printf("%15.12f %15.12f\n",(ko[spin][i]+0.15)*27.2113845,sum);
    }
  }
  
}  // the end of 'main' 




void Draw_Gcube(double coef[6])
{
  int    i,j,k;       // variable for loop
  double x,y,z;       // coordinate 
  int  atomnum = 7;   // the number of atoms
  double X=10.;       /*  each length of the box which contain the whole wave ftn.
		          Unit : Angstrom    */
  double Y=10.;
  double Z=10.;
  double Origin_x = 0. ; // X/2 ;
  double Origin_y = 0. ; // Y/2 ;
  double Origin_z = 0. ; // Z/2 ;

  int Nx=101;            // the number of grid(values) in each direction
  int Ny=101;
  int Nz=101; 

  double dx=X/(Nx-1.);     // mesh of each direction
  double dy=Y/(Ny-1.);
  double dz=Z/(Nz-1.);
  


  /*** variables for constructing the wave ftn. ***/ 
  double a;
  double r;
  double R32;
  double d1,d2,d3,d4,d5;
  double Psi1;    // Psi2, Psi3, Psi4, Psi5;


  FILE *fout;

  fout = fopen("./Mn12_.cube", "w");

  fprintf(fout,"Title of ..\n");
  fprintf(fout,"Title of ..\n");
  fprintf(fout,"%5d   %+17.13E  %+17.13E  %+17.13E\n", atomnum,Origin_x,Origin_y,Origin_z);
  
  fprintf(fout, "%5d   %+17.13E  %+17.13E  %+17.13E\n", Nx, dx, 0., 0.);
  fprintf(fout, "%5d   %+17.13E  %+17.13E  %+17.13E\n", Ny, 0., dy, 0.);
  fprintf(fout, "%5d   %+17.13E  %+17.13E  %+17.13E\n", Nz, 0., 0., dz);


  /* Mn : atomic number, core charge number, x, y, z - coordinate */

  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 25, 2., X/2., Y/2., Z/2.);

  /* O : atomic number, core charge number, x, y, z - coordinate */

  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2.,
           X/2.-(Gxyz[Choo_atom[0]][1]-Gxyz[Choo_atom[1]][1]),
           Y/2.-(Gxyz[Choo_atom[0]][2]-Gxyz[Choo_atom[1]][2]),
           Z/2.-(Gxyz[Choo_atom[0]][3]-Gxyz[Choo_atom[1]][3]));

  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2.,
           X/2.-(Gxyz[Choo_atom[0]][1]-Gxyz[Choo_atom[2]][1]),
           Y/2.-(Gxyz[Choo_atom[0]][2]-Gxyz[Choo_atom[2]][2]),
           Z/2.-(Gxyz[Choo_atom[0]][3]-Gxyz[Choo_atom[2]][3]));

  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2.,
           X/2.-(Gxyz[Choo_atom[0]][1]-Gxyz[Choo_atom[3]][1]),
           Y/2.-(Gxyz[Choo_atom[0]][2]-Gxyz[Choo_atom[3]][2]),
           Z/2.-(Gxyz[Choo_atom[0]][3]-Gxyz[Choo_atom[3]][3]));

  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2.,
           X/2.-(Gxyz[Choo_atom[0]][1]-Gxyz[Choo_atom[4]][1]),
           Y/2.-(Gxyz[Choo_atom[0]][2]-Gxyz[Choo_atom[4]][2]),
           Z/2.-(Gxyz[Choo_atom[0]][3]-Gxyz[Choo_atom[4]][3]));

  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2.,
           X/2.-(Gxyz[Choo_atom[0]][1]-Gxyz[Choo_atom[5]][1]),
           Y/2.-(Gxyz[Choo_atom[0]][2]-Gxyz[Choo_atom[5]][2]),
           Z/2.-(Gxyz[Choo_atom[0]][3]-Gxyz[Choo_atom[5]][3]));

  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2.,
           X/2.-(Gxyz[Choo_atom[0]][1]-Gxyz[Choo_atom[6]][1]),
           Y/2.-(Gxyz[Choo_atom[0]][2]-Gxyz[Choo_atom[6]][2]),
           Z/2.-(Gxyz[Choo_atom[0]][3]-Gxyz[Choo_atom[6]][3]));

  /* O17 */
  /*
  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2., X/2.-0.119501, Y/2.+1.577760, Z/2.-1.039353);
  */

  /* O22 */
  /*
  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2., X/2.+1.775198, Y/2.+0.381018, Z/2.+0.563654);
  */

  /* O37 */
  /*
  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2., X/2.-0.846899, Y/2.+0.869413, Z/2.+1.754141);
  */

  /* O45 */
  /*
  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2., X/2.+0.041566, Y/2.-1.693799, Z/2.+1.020771);
  */

  /* O53 */
  /*
  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2., X/2.-1.844473, Y/2.-0.602702, Z/2.-0.526490);
  */

  /* O57 */
  /*
  fprintf(fout,"%3d  %+15.12E   %+15.12E  %+15.12E  %+15.12E\n"
	  , 8, 2., X/2.+0.609629, Y/2.-1.179424, Z/2.-1.726887);
  */

  for(k=0; k<Nz; k++)
    {
      for(j=0; j<Ny; j++)
	{
	  for(i=0; i<Nx; i++)
	    {
	      
	      /*** construction of the wave ftn - Psi1 ***/

	      x= dx*i - X/2 ;
	      y= dy*j - Y/2 ;
	      z= dz*k - Z/2 ;

	      a = sqrt( 15 / (4*3.141592654) );  // nomalization const.
	      r = sqrt( x*x + y*y + z*z );
	      R32 =  exp(-r/1.);   // R(r)-ftn  

	      d1 = (3*z*z-r*r)/(2*1.73205081);
	      d2 = (x*x-y*y)/(2.);
	      d3 = (x*y);
	      d4 = (x*z);
	      d5 = (y*z);
	      
	      Psi1 = a*R32*(coef[1]*d1 + coef[2]*d2 + coef[3]*d3 + coef[4]*d4 + coef[5]*d5);

	      fprintf(fout, "%+13.6e  ", Psi1);

	    }
	  fprintf(fout,"\n");
	}
    }

  fclose(fout);

}




void Eigen_lapack(double **a, double *ko, int n0)
{
 /* input:  n;
    input:  a[n][n];  matrix A

    output: a[n][n]; eigevectors
    output: ko[n];   eigenvalues  */
    
  static char *name="Eigen_lapack";

  char  *JOBZ="V";
  char  *RANGE="A";
  char  *UPLO="L";

  INTEGER n=n0;
  INTEGER LDA=n;
  double VL,VU; /* dummy */
  INTEGER IL,IU; /* dummy */
  double ABSTOL=1.0e-10;
  INTEGER M;

  double *A,*Z;
  INTEGER LDZ=n;
  INTEGER LWORK;
  double *WORK;
  INTEGER *IWORK;

  INTEGER *IFAIL, INFO;

  INTEGER i,j;

  A=(double*)malloc(sizeof(double)*n*n);
  Z=(double*)malloc(sizeof(double)*n*n);

  LWORK=n*8;
  WORK=(double*)malloc(sizeof(double)*LWORK);
  IWORK=(INTEGER*)malloc(sizeof(INTEGER)*n*5);
  IFAIL=(INTEGER*)malloc(sizeof(INTEGER)*n);
 

  for (i=0;i<n;i++) {
    for (j=0;j<n;j++) {
       A[i*n+j]= a[i+1][j+1];
    }
  }

#if 0
  printf("A=\n");
  for (i=0;i<n;i++) {
    for (j=0;j<n;j++) {
       printf("%f ",A[i*n+j]);
    }
    printf("\n");
  }
  fflush(stdout);
#endif

  F77_NAME(dsyevx,DSYEVX)( JOBZ, RANGE, UPLO, &n, A, &LDA, &VL, &VU, &IL, &IU,
           &ABSTOL, &M, ko, Z, &LDZ, WORK, &LWORK, IWORK,
           IFAIL, &INFO ); 

  /* store eigenvectors */
  for (i=0;i<n;i++) {
    for (j=0;j<n;j++) {
     /*  a[i+1][j+1]= Z[i*n+j]; */
       a[j+1][i+1]= Z[i*n+j];

    }
  }

  /* shift ko by 1 */
  for (i=n;i>=1;i--){
    ko[i]= ko[i-1];
  }

  if (INFO>0) {
     printf("\n%s: error in dsyevx_, info=%d\n\n",name,INFO);
  }
  if (INFO<0) {
     printf("%s: info=%d\n",name,INFO);
     exit(10);
  }
   
  free(IFAIL); free(IWORK); free(WORK); free(Z); free(A);

}


void read_coordinates(char *argv[])
{
  static FILE *fp;
  static char Species[20];
  static double tmp0,tmp1;
  static int i,j;

  Gxyz = (double**)malloc(sizeof(double*)*(atomnum+1)); 
  for (i=0; i<=atomnum; i++){
    Gxyz[i] = (double*)malloc(sizeof(double)*60); 
  }

  if ((fp = fopen(argv[2],"r")) != NULL){
    printf("Read the coordinates file (%s)\n",argv[2]);

    /****************************************************
                         open the file
    ****************************************************/

    input_open(argv[2]);

    /****************************************************
                         read data
    ****************************************************/

    if (fp=input_find("<Atoms.SpeciesAndCoordinates") ) {
      for (i=1; i<=atomnum; i++){
        fscanf(fp,"%i %s %lf %lf %lf %lf %lf",&j,&Species,
	       &Gxyz[i][1],&Gxyz[i][2],&Gxyz[i][3],&tmp0,&tmp1);

        if (i!=j){
          printf("Format error of the sequential number %i in <Atoms.SpeciesAndCoordinates\n",j);
          exit(0);
        }
      }

      if (!input_last("Atoms.SpeciesAndCoordinates>")) {
        /* format error */
        printf("Format error for Atoms.SpeciesAndCoordinates\n");
        exit(0);
      }
    }

    /****************************************************
                       input_close
    ****************************************************/

    input_close();
  }
  else {
    printf("Failure of reading the coordinates file (%s).\n",argv[2]);
  }

  for (i=1; i<=atomnum; i++){
    printf("i=%i  coodinates=%15.12f %15.12f %15.12f\n",
	   i,Gxyz[i][1],Gxyz[i][2],Gxyz[i][3]);
  }

}