#include<iostream>
#include<fstream>
#include<math.h>
#include<iomanip>

//CLASS READS ATOMIC POTENTIALS FROM ATPOT_C AND ATPOT_O DIRECTORIES
//ATOMIC POTENTIALS ARE OBTAINED FROM UNRESTRICTED CALCULATIONS USING A (2x2) CELL 

/****************************************************************************************************************************/
/* PROGRAM BY GERNOT FÜCHSEL, APRIL 2013                                                                                    */
/* UNIVERSITY OF POTSDAM                                                                                                    */
/* THEORETICAL CHEMISTRY GROUP PETER SAALFRANK                                                                              */
/*                                                                                                                          */
/* MODIFIED JANUARY 2016 AT LEIDEN UNIVERSITY GROPU OF GEERT-JAN KROES                                                      */
/* FOR STEPPED SURFAACE Cu(211)                                                                                             */
/*                                                                                                                          */
/* FOR SURFACE SCATTERING PROCESSES OF                                                                                      */
/* HOMONUCLEAR MOLECULES AT SURFACES WITH CS SYMMETRY AT TOP POSITION                                                       */
/* THE CLASS ATOMIC.h INTERPOLATES THE POTENTIAL OF AN ATOM INCORPORATING CS SYMMETRY                                       */
/****************************************************************************************************************************/

using namespace std;


//extern void invert_matrix(double* MATRIX, int M);
extern  void spline_init(double*, double*, int, double*);
extern  void spline_interpolation(double*, double*, double*, int, double, double*, double* );

//LAPACK ROUTINES
extern "C"{void dgetri_(int* M, double* A, int* LDA, int* IPIV, double* WORK, int* LWORK, int* INFO);}
extern "C"{void dgetrf_(int* M, int* N, double* A, int* LDA, int* IPIV, int* INFO );}


class ATOMIC{
 private:
  const static double LX = 6.955149254;   // lattice constant along X
  const static double LY = 2.839427793;   // lattice constant along X
  const static int nvrepH       = 163;             // Number of vrep potential points to be read
  const static int NH           = 109;             // Number of vrep potential points to be read
  int n_clust    ;
  double offset  ;                                //asymtotic potential value (value at large Z)
  double* cluster_pos;                            // coordinates of cluster atoms

  const static int Nbasis = 27; //number of basis functions to be used in the u, v interpolation
  double* INVERSE_Matrix;
  double* vpot_site; double* dvd_site;
  double* CELL_CLUSTER;
  int NLAYER1; int NLAYER2; int NLAYER3;
  double* LAYER1; double* LAYER2; double* LAYER3;



  //Construct Cluster for reference functions
  void BUILD_CLUSTER();
  // calculates the reference potential functions for corrugation reduced scheme and its derivates
  // V(x,y,z) = I(x,y,z) + sum_i Q_i(x,y,z)
  void REF_POT(double x, double y, double z, double* v, double* dvdx, double* dvdy, double* dvdz, int mode);
  void read_files(ifstream &inp, char* file, double* V, double* moms, double* Z, double X, double Y, double offset); 
  //Return to ORIGIN
  void SET_TO_ORIGIN(double* x, double* y); // returns position back to unit cell
  void CALC_ATOMIC_UV(double u0, double v0, double* vpot_site, double* dvd_site, double* vpot, double* dvdu, double* dvdv, double* dvdz);
  //INITIALIZE ARRAYS
  void INIT_ARRAYS();
  double* vrepH   ; // repulsive potential at (u,v) = (0,0); C_ATOM
  double* rvrepH  ; // distance to top position

  double* rvrepH2; double* vrepH2; 
  const static int nvrepH2 = 163;//21;
  double* moms_vrepH2;
  
  //arrays for spline interpolation
  double* moms_vrepH; 
  double* ZH;
  //POTENTIAL VALUES AND SPLINE COEFFICIENTS ARE STORED IN FULL_ARRAY
  double* FULL_ARRAY; double* BASIS_FUNC; double* INVERSE;
  double* DBASIS_DU; double* DBASIS_DV;
  double gux; double guy;
  void set_basis_val(double u, double v);
  double get_basis(int i){return(BASIS_FUNC[i]);};
  double get_dbasis_du(int i){return(DBASIS_DU[i]);};
  double get_dbasis_dv(int i){return(DBASIS_DV[i]);};
  void invert_matrix(double* MATRIX, int M);
  

 public:
  void INIT_ATOMIC(char* direc); 
  void CALC_ATOMIC(double* pos, double* dvd, double* v1);
}; // end member definition
void ATOMIC::invert_matrix(double* MATRIX, int M){
   //FIRST DETERMINE THE PIVOT ELEMENTS
   int INFO; int LDA = M; int* IPIV = new int[M];
   int N = M;
   int LWORK = M;
   double* WORK = new double[LWORK];

   dgetrf_(&M, &N, MATRIX, &LDA, IPIV, &INFO );
   if(INFO != 0) std::cerr << "ERROR: PIVOT- ELEMENTS COULD NOT BE DETERMINED   INFO = " << INFO << "\n";

   //INVERT THE MATRIX
   dgetri_(&M, MATRIX, &LDA, IPIV, WORK, &LWORK, &INFO);
   if(INFO != 0) std::cerr << "ERROR: MATRIX COULD NOT BE INVERTED \n";

   delete IPIV; delete WORK;
};
void  ATOMIC::set_basis_val(double u, double v){
//EVALUATES THE BASIS FUNCTIONS
  BASIS_FUNC[0]  = 1.;                                   
  BASIS_FUNC[1]  = cos(u);
  BASIS_FUNC[2]  = cos(2.*u);
  BASIS_FUNC[3]  = cos(v);
  BASIS_FUNC[4]  = cos(u+v) + cos(u-v);
  BASIS_FUNC[5]  = cos(3.*u);
  BASIS_FUNC[6]  = cos(2.*u+v) + cos(2.*u-v);
  BASIS_FUNC[7]  = cos(3.*u+v) + cos(3.*u-v);
  BASIS_FUNC[8]  = cos(4.*u);
  BASIS_FUNC[9]  = cos(4.*u+v) + cos(4.*u-v);
  BASIS_FUNC[10] = cos(2.*v);
  BASIS_FUNC[11] = cos(u+2*v) + cos(u-2.*v);
  BASIS_FUNC[12] = cos(2.*(u+v)) + cos(2.*(u-v));
  BASIS_FUNC[13] = cos(3.*u+2.*v) + cos(3.*u-2.*v);
  BASIS_FUNC[14] = cos(4.*u + 2.*v) + cos(4.*u-2.*v);
  BASIS_FUNC[15] = sin(u);
  BASIS_FUNC[16] = sin(2.*u);
  BASIS_FUNC[17] = sin(u+v) + sin(u-v);
  BASIS_FUNC[18] = sin(3.*u);
  BASIS_FUNC[19] = sin(2.*u+v) + sin(2.*u-v);
  BASIS_FUNC[20] = sin(3.*u+v) + sin(3.*u-v);
  BASIS_FUNC[21] = sin(4.*u);
  BASIS_FUNC[22] = sin(4.*u+v) + sin(4.*u-v);
  BASIS_FUNC[23] = sin(u+2.*v) + sin(u-2*v);
  BASIS_FUNC[24] = sin(2.*(u+v)) + sin(2.*(u-v));
  BASIS_FUNC[25] = sin(3.*u + 2.*v) + sin(3.*u - 2.*v);
  BASIS_FUNC[26] = sin(4.*u+2.*v) + sin(4.*u-2.*v);

//ANALYTICAL DERIVATIVES CAN BE IMPLEMENTED AT A LATER POINT
  
  DBASIS_DU[0]  = 0.;                     DBASIS_DV[0]  = 0.;
  DBASIS_DU[1]  = -BASIS_FUNC[15];        DBASIS_DV[1]  = 0.;
  DBASIS_DU[2]  = -2.*BASIS_FUNC[16];     DBASIS_DV[2]  = 0.;
  DBASIS_DU[3]  = 0.;                     DBASIS_DV[3]  = -sin(v);
  DBASIS_DU[4]  = -BASIS_FUNC[17];        DBASIS_DV[4]  = -sin(u+v) + sin(u-v);
  DBASIS_DU[5]  = -3.*BASIS_FUNC[18];     DBASIS_DV[5]  = 0.;
  DBASIS_DU[6]  = -2.*BASIS_FUNC[19];     DBASIS_DV[6]  = -sin(2.*u+v) + sin(2.*u-v);
  DBASIS_DU[7]  = -3.*BASIS_FUNC[20];     DBASIS_DV[7]  = -sin(3.*u+v) + sin(3.*u-v);
  DBASIS_DU[8]  = -4.*BASIS_FUNC[21];     DBASIS_DV[8]  = 0.;
  DBASIS_DU[9]  = -4.*BASIS_FUNC[22];     DBASIS_DV[9]  = -sin(4.*u+v) + sin(4.*u-v);
  DBASIS_DU[10] = 0.;                     DBASIS_DV[10] = -2.*sin(2.*v);
  DBASIS_DU[11] =    -BASIS_FUNC[23];     DBASIS_DV[11] = 2.*(-sin(u+2*v)     + sin(u-2.*v)  );
  DBASIS_DU[12] = -2.*BASIS_FUNC[24];     DBASIS_DV[12] = 2.*(-sin(2.*(u+v))  + sin(2.*(u-v)));
  DBASIS_DU[13] = -3.*BASIS_FUNC[25];     DBASIS_DV[13] = 2.*(-sin(3.*u+2.*v) + sin(3.*u-2.*v));
  DBASIS_DU[14] = -4.*BASIS_FUNC[26];     DBASIS_DV[14] = 2.*(-sin(4.*u+2.*v) + sin(4.*u-2.*v));
  DBASIS_DU[15] =    BASIS_FUNC[1];       DBASIS_DV[15] = 0.;
  DBASIS_DU[16] = 2.*BASIS_FUNC[2];       DBASIS_DV[16] = 0.;        
  DBASIS_DU[17] =    BASIS_FUNC[4];       DBASIS_DV[17] = cos(u+v) - cos(u-v);
  DBASIS_DU[18] = 3.*BASIS_FUNC[5];       DBASIS_DV[18] = 0.;
  DBASIS_DU[19] = 2.*BASIS_FUNC[6];       DBASIS_DV[19] = cos(2.*u+v) - cos(2.*u-v);
  DBASIS_DU[20] = 3.*BASIS_FUNC[7];       DBASIS_DV[20] = cos(3.*u+v) - cos(3.*u-v);
  DBASIS_DU[21] = 4.*BASIS_FUNC[8];       DBASIS_DV[21] = 0.;
  DBASIS_DU[22] = 4.*BASIS_FUNC[9];       DBASIS_DV[22] = cos(4.*u+v) - cos(4.*u-v);        
  DBASIS_DU[23] =    BASIS_FUNC[11];      DBASIS_DV[23] = 2.*(cos(u+2.*v) - cos(u-2*v));
  DBASIS_DU[24] = 2.*BASIS_FUNC[12];      DBASIS_DV[24] = 2.*(cos(2.*(u+v))  - cos(2.*(u-v)));
  DBASIS_DU[25] = 3.*BASIS_FUNC[13];      DBASIS_DV[25] = 2.*(cos(3.*u+2.*v) - cos(3.*u-2.*v));
  DBASIS_DU[26] = 4.*BASIS_FUNC[14];      DBASIS_DV[26] = 2.*(cos(4.*u+2.*v) - cos(4.*u-2.*v));
     

};

/**************************************************************************************************/
void ATOMIC::SET_TO_ORIGIN(double* x, double* y){
  double u0 = (*x)/LX;
  double v0 = (*y)/LY; 
  if(u0 != fabs(u0)) u0 += 1;
  if(v0 != fabs(v0)) v0 += 1;
  u0 -= (int)u0;
  v0 -= (int)v0;
  *x = u0*LX;
  *y = v0*LY;
};
/**************************************************************************************************/

void ATOMIC::BUILD_CLUSTER(){
//defines the bulk atom positions that are needed to calculate
//the reference potentials in the corrugation reduced procedure

  int nu = 2; int nv = 3; //number of shift vectors from -u to u, and from -v to v
  int num_atoms = 8;     //number of atoms in a 1x1 cell
  n_clust     = num_atoms*(2*nu+1)*(2*nv+2);   //total number of atoms in the cluster

  NLAYER1 = 3*(2*nu+1)*(2*nv+2);
  NLAYER2 = 3*(2*nu+1)*(2*nv+2);
  NLAYER3 = 2*(2*nu+1)*(2*nv+2);

  LAYER1 = new double[3*NLAYER1];
  LAYER2 = new double[3*NLAYER2];
  LAYER3 = new double[3*NLAYER3];

  CELL_CLUSTER = new double[num_atoms*3];
  //ATOM 1 is fixed at the origin of the COORDINATE SYSTEM
  //CARTESIAN COORDINATES (DIRECT COORDINATES IN XY MULTIPLIED WITH LX OR LY)
  //FIRST LAYER
  CELL_CLUSTER[0]  = 0.0            ; CELL_CLUSTER[1]  = 0.0    ;  CELL_CLUSTER[2]  = 0.0;
  CELL_CLUSTER[3]  = 0.6745912857*LX; CELL_CLUSTER[4]  = 0.5*LY ;  CELL_CLUSTER[5]  = -0.60063417815885;
  CELL_CLUSTER[6]  = 0.3540812621*LX; CELL_CLUSTER[7]  = 0.0    ;  CELL_CLUSTER[8]  = -1.26906089571607;
  //SECOND LAYER
  CELL_CLUSTER[9]  = 0.0313933192*LX; CELL_CLUSTER[10] = 0.5*LY ;  CELL_CLUSTER[11] = -2.32097067828095;
  CELL_CLUSTER[12] = 0.6892865112*LX; CELL_CLUSTER[13] = 0.0    ;  CELL_CLUSTER[14] = -3.05526814621046;
  CELL_CLUSTER[15] = 0.3637807448*LX; CELL_CLUSTER[16] = 0.5*LY ;  CELL_CLUSTER[17] = -3.88075461289106;
  //THIRD LAYER
  CELL_CLUSTER[18] = 0.0289107484*LX; CELL_CLUSTER[19] = 0.0    ;  CELL_CLUSTER[20] = -4.67980858324286;
  CELL_CLUSTER[21] = 0.6936261601*LX; CELL_CLUSTER[22] = 0.5*LY ;  CELL_CLUSTER[23] = -5.49232563522636;
  double x,y,z;
  int count = 0;
  for(int i = -nu; i <= nu; i++){
    for(int j = -nv; j <= nv+1; j++){
      for(int k = 0; k < 3; k++){ //first layer
        x = CELL_CLUSTER[3*k]   + i*LX;
        y = CELL_CLUSTER[3*k+1] + j*LY;
        z = CELL_CLUSTER[3*k+2]; 
        LAYER1[(count + k)*3 + 0] = x;
        LAYER1[(count + k)*3 + 1] = y;
        LAYER1[(count + k)*3 + 2] = z;
//        cout << "HERE : " << k << "\n";
      }
      for(int k = 0; k < 3; k++){ //second Layer layer
        x = CELL_CLUSTER[3*(k+3)]   + i*LX;
        y = CELL_CLUSTER[3*(k+3)+1] + j*LY;
        z = CELL_CLUSTER[3*(k+3)+2];
        LAYER2[(count + k)*3 + 0] = x;
        LAYER2[(count + k)*3 + 1] = y;
        LAYER2[(count + k)*3 + 2] = z;
      }
      for(int k = 0; k < 2; k++){ //THIRD LAYER
        x = CELL_CLUSTER[3*(k+6)]   + i*LX;
        y = CELL_CLUSTER[3*(k+6)+1] + j*LY;
        z = CELL_CLUSTER[3*(k+6)+2];
        LAYER3[count*2 + k*3 + 0] = x;
        LAYER3[count*2 + k*3 + 1] = y;
        LAYER3[count*2 + k*3 + 2] = z;
      }
      count +=3;
    }
  }
  
//WRTE CLUSTER OUT

ofstream out; out.open("cluster.xyz");
out << n_clust << "\n\n";
for(int i = 0; i < NLAYER1; i++){
  out << "Pt " << setprecision(13)<< showpoint << LAYER1[i*3] << " " << LAYER1[i*3 + 1] << " " << LAYER1[i*3 + 2] << "\n";
}
for(int i = 0; i < NLAYER2; i++){
  out << "Cu " << LAYER2[i*3] << " " << LAYER2[i*3 + 1] << " " << LAYER2[i*3 + 2] << "\n";
}
for(int i = 0; i < NLAYER3; i++){
  out << "Ru " << LAYER3[i*3] << " " << LAYER3[i*3 + 1] << " " << LAYER3[i*3 + 2] << "\n";
}

};

/**************************************************************************************************/
void ATOMIC::REF_POT(double x, double y, double z, double* v, double* dvdx, double* dvdy, double* dvdz, int mode){

  //REF_POT calculates the reference potentials und the corresponding derivates 
  //needed in the corrugation reduced procedure

  double r       = 0.; double vpot    = 0.; 
  double dvpotdx = 0.; double dvpotdy = 0.; double dvpotdz = 0.;
  double svpot, dvd  ; double X,Y,Z;

  if(mode == 0){ 

  for(int i = 0; i < NLAYER1; i++){
    X = (x - LAYER1[i*3])  ;
    Y = (y - LAYER1[i*3+1]);
    Z = (z - LAYER1[i*3+2]);
    r = sqrt(X*X+Y*Y+Z*Z);
    if(r<=rvrepH[nvrepH-1]){
      spline_interpolation(rvrepH, vrepH, moms_vrepH, nvrepH, r, &svpot, &dvd);
      vpot    += svpot;
      dvpotdx += X*dvd/r;
      dvpotdy += Y*dvd/r;
      dvpotdz += Z*dvd/r;
    }
  }//end for  loop
  }
  else if(mode == 1){
  for(int i = 0; i < NLAYER2; i++){
    X = (x - LAYER2[i*3])  ;
    Y = (y - LAYER2[i*3+1]);
    Z = (z - LAYER2[i*3+2]);
    r = sqrt(X*X+Y*Y+Z*Z);
    if(r<=rvrepH2[nvrepH2-1]){
      spline_interpolation(rvrepH2, vrepH2, moms_vrepH2, nvrepH2, r, &svpot, &dvd);
      vpot    += svpot;
      dvpotdx += X*dvd/r;
      dvpotdy += Y*dvd/r;
      dvpotdz += Z*dvd/r;
    }
  }//end for  loop
  }
  else if(mode == 2){
  for(int i = 0; i < NLAYER3; i++){
    X = (x - LAYER3[i*3])  ;
    Y = (y - LAYER3[i*3+1]);
    Z = (z - LAYER3[i*3+2]);
    r = sqrt(X*X+Y*Y+Z*Z);
    if(r<=rvrepH2[nvrepH2-1]){
      spline_interpolation(rvrepH2, vrepH2, moms_vrepH2, nvrepH2, r, &svpot, &dvd);
      vpot    += svpot;
      dvpotdx += X*dvd/r;
      dvpotdy += Y*dvd/r;
      dvpotdz += Z*dvd/r;
    }
  }//end for  loop
  }
  else{ cerr << "ERROR IN REFERENCE POT. PLEASE PROVIDE A MODE\n";}

  *v    = vpot;
  *dvdx = dvpotdx;
  *dvdy = dvpotdy;
  *dvdz = dvpotdz;


};


/**************************************************************************************************/
void ATOMIC::INIT_ARRAYS(){
  //IN THIS ROUTINE THE MOST IMPORTANT CONSTANTS AND ARRAYS ARE INITIALIZED. 
  //FURTHER, A SET OF LINEAR EQUATIONS IS SOLVED TO OBTAIN THE FOURIER COEFFICIENTS
  //FOR THE INTERPOLATION ALONG U AND V.
  //THE COEFFICIENTS ARE STORED IN "INVERSE_Matrix"
  

  //INITIALIZE SOME CONSTANTS OFTEN USED IN THE CLASS 
  //AND SOME ARRAYS
  vrepH         = new double[nvrepH]; // repulsive potential at (u,v) = (0,0); C_ATOM
  rvrepH        = new double[nvrepH]; // distance to top position
  
  //arrays for spline interpolation
  moms_vrepH    = new double[nvrepH];
  ZH = new double[NH];

  rvrepH2 = new double[nvrepH]; 
  vrepH2 = new double[nvrepH]; 
  moms_vrepH2 = new double[nvrepH];

  FULL_ARRAY = new double[Nbasis*NH*2]; //store the potential values and the spline coefficients 
  BASIS_FUNC = new double[Nbasis];
  DBASIS_DU  = new double[Nbasis];
  DBASIS_DV  = new double[Nbasis];
  INVERSE_Matrix = new double[Nbasis*Nbasis];

  vpot_site = new double[Nbasis]; 
  dvd_site  = new double[Nbasis];

  gux = 2.*M_PI/LX;
  guy = 2.*M_PI/LY;

  //SETUP THE INVERSE_MATRIX
  double u , v;
  //DEFINITION OF THE GRID-POINTS
  for(int i = 0; i < 3; i++){
    v = 0.25*i * LY * guy;
    for(int j = 0; j < 9; j++){
      u = j * (1./9. * LX * gux);  //this respects that the positions are defined equidistantly on the skewed stepped plane
      set_basis_val(u,v); 
      for(int k = 0; k < Nbasis; k++){
        double val = get_basis(k);
        if(fabs(val) < 1e-12) val = 0.; //for numerics
        INVERSE_Matrix[(i*9+j)*Nbasis + k] = val; //store according to: rows are positions, columns are functions
      }
    }
  }  
  invert_matrix(INVERSE_Matrix, Nbasis);
/*
  for(int i = 0; i < Nbasis; i++){
    for(int j = 0; j < Nbasis; j++){ cout << setprecision(7) << fixed << setw(10) << INVERSE_Matrix[i*Nbasis+j] << " ";
    }
    cout << "\n";
  }  
*/ 
  
  

};
/**************************************************************************************************/
void ATOMIC::INIT_ATOMIC(char* direc){
  char file_name[1024];
  sprintf(file_name,"%s/vrep.dat",direc); 
  ifstream ivrep; ivrep.open(file_name);

  cout << file_name << "\n";

  //INITIALZE THE ARRAYS
  INIT_ARRAYS();
  //CONSTRUCT THE REFERENCE CLUSTER FOR REFERENCE POTENTIALS
  BUILD_CLUSTER();
  // read repulsive potential of the C-ATOM
  for(int i = 0; i < nvrepH; i++){ 
    ivrep >> rvrepH[i] >> vrepH[i];// >>  muell;
    rvrepH2[i] = rvrepH[i];
  }
  offset = vrepH[nvrepH-1];
  //make sure that for r -> infinity vrep_C = 0.;
  for(int i = 0; i < nvrepH; i++)  vrepH[i] -= offset;
  for(int i = 0; i < nvrepH; i++)  vrepH2[i] = vrepH[i];

  spline_init(rvrepH, vrepH, nvrepH, moms_vrepH);
  spline_init(rvrepH2, vrepH2, nvrepH2, moms_vrepH2);  

  ivrep.close();

  double X, Y;
  ifstream inpat  ;
  //WORK STORES THE POTENTIAL AND WORK2 STORES THE SPLINE COEFFICIENTS
  double* work = new double[NH]; double* work2 = new double[NH];
  int count = 0; 

/*
  double dum1, dum2;//, dvdx,dvdy, dvdz, vref;
  //INITIALZE THE NEXT REFERENCE POTENTIAL
  sprintf(file_name,"%s/vrep_top3.dat",direc);
  cout << "# " << file_name << "\n";
  ivrep.open(file_name);
  X = 3./9.*LX; Y = 0.;
  for(int i = 0; i < nvrepH2; i++){
    ivrep >> dum1 >> dum2;
    rvrepH2[i] = (dum1 + 2.87994); 
//    REF_POT(X, Y, dum1, &vref, &dvdx, &dvdy, &dvdz,0);
    vrepH2[i]  = dum2;   
//    cout << dum1 << " " <<  rvrepH2[i] << " " << vrepH2[i] << "\n";
  }
  for(int i = 0; i < nvrepH2; i++) vrepH2[i] -= vrepH2[nvrepH2-1];
  //make sure that for r -> infinity vrep_C = 0.;
  spline_init(rvrepH2, vrepH2, nvrepH2, moms_vrepH2);
  ivrep.close();
*/

  
  double muell;
  //Z-VALUES FOR THE 1D POTENTIAL CUTS
  sprintf(file_name,"%s/TOP_0.dat",direc);
  inpat.open(file_name);
  for(int i = 0; i < NH; i++){inpat >>  ZH[i] >> muell;}
  inpat.close();
  //FIRST READ THE TOP LINE
  for(int i = 0; i < 9; i++){
    sprintf(file_name,"%s/TOP_%d.dat",direc, i);               
    inpat.open(file_name);
    X = (LX/9.)*i;      Y = 0.;
    read_files(inpat, file_name, work, work2, ZH, X, Y, offset);
    inpat.close();
    for(int j = 0; j < NH; j++) FULL_ARRAY[count + j] = work[j];
    count += NH;
    for(int j = 0; j < NH; j++) FULL_ARRAY[count + j] = work2[j];
    count += NH;
  }
 
  //SECOND READ THE T2B LINE
  for(int i = 0; i < 9; i++){
    sprintf(file_name,"%s/T2B_%d.dat",direc, i);
    inpat.open(file_name);
    X = (LX/9.)*i;      Y = LY/4.;
    read_files(inpat, file_name, work, work2, ZH, X, Y, offset);
    inpat.close();
    for(int j = 0; j < NH; j++) FULL_ARRAY[count + j] = work[j];
    count += NH;
    for(int j = 0; j < NH; j++) FULL_ARRAY[count + j] = work2[j];
    count += NH;
  }
  
  //THIRD READ THE BRG LINE
  for(int i = 0; i < 9; i++){
    sprintf(file_name,"%s/BRG_%d.dat",direc, i);
    inpat.open(file_name);
    X = (LX/9.)*i;      Y = 2.*LY/4.;
    read_files(inpat, file_name, work, work2, ZH, X, Y, offset);
    inpat.close();
    for(int j = 0; j < NH; j++) FULL_ARRAY[count + j] = work[j];
    count += NH;
    for(int j = 0; j < NH; j++) FULL_ARRAY[count + j] = work2[j];
    count += NH;
  }
  delete work; delete work2; 

  ofstream out; out.open("top_crp.info");
  for(int i = 0; i < NH; i++){
    out << ZH[i] << " "; 
    for(int j = 0; j < 9; j++) out << setw(14) << FULL_ARRAY[2*NH*j + i] << "  ";
    out << "\n";
  }
  out.close();
  out.open("t2b_crp.info");
  for(int i = 0; i < NH; i++){
    out << ZH[i] << " ";
    for(int j = 9; j < 18; j++) out << setw(14) << FULL_ARRAY[2*NH*j + i] << "  ";
    out << "\n";
  }
  out.close();
  out.open("brg_crp.info");
  for(int i = 0; i < NH; i++){
    out << ZH[i] << " ";
    for(int j = 18; j < 27; j++) out << setw(14) << FULL_ARRAY[2*NH*j + i] << "  ";
    out << "\n";
  }
  out.close();

//  cout << "INIT DONE \n";

 
}; // definition of INIT_POT_C
void ATOMIC::read_files(ifstream &inp, char* file, double* V, double* moms, double* Z, double X, double Y, double offset){
  double muell; double vref; double dvdx, dvdy, dvdz; double vref2 = 0.;
  //CHECK STATUS OF THE FILE
  if ( inp.rdstate() != 0) cerr << "FATAL ERROR IN FILE " << file << "\n";
  for(int i = 0; i < NH; i++){
    inp >> muell >> V[i]; 
    REF_POT(X, Y, Z[i], &vref, &dvdx, &dvdy, &dvdz, 0);
    vref2 = vref;
//    cout << muell <<  " " << vref  << " " ;
    REF_POT(X, Y, Z[i], &vref, &dvdx, &dvdy, &dvdz, 1);
    vref2 += vref;
//    cout << vref << "\n" ;
    REF_POT(X, Y, Z[i], &vref, &dvdx, &dvdy, &dvdz, 2);
    vref2 += vref;
//    cout << vref << "\n";
    V[i] -= vref2+offset;
//   V[i] -= offset; 
  }
  spline_init(Z, V, NH, moms);
};

/**************************************************************************************************/
void ATOMIC::CALC_ATOMIC(double* pos_cart, double* dvd_cart, double* vb1){
  //3D POTENTIAL ROUTINE (u;v;z)
  //TRIGONOMETRIC INTERPOLATION ALONG X,Y.
  //SPLINE INTERPOLATION ALONG Z FOR 1D REFERENCE FUNCTIONS
  //SINCE X,Y ARE NOT THE COORDINATES TO WORK WITH,
  //THE COORDINATE SYSTEM IS CHANGED TO DIRECT LATTICE VECTORS U,V. 
  double vref = 0.; double vpot = 0.; double dvd = 0.;    
  double dvdx1, dvdy1, dvdz11; 
  double dvdu1 = 0.; double dvdv1 = 0.; double dvdz1 = 0.;
  double x1 = pos_cart[0]*0.52917721; 
  double y1 = pos_cart[1]*0.52917721; 
  double z1 = pos_cart[2]*0.52917721;


  for(int i = 0; i < Nbasis; i++) vpot_site[i] = dvd_site[i] = 0.;
  double u0 = x1*gux; double v0 = y1*guy;
  //need two pointers
  double* pot; double* moms; int count = 0;
  if(z1 <= ZH[NH-1]){
    //NOTE IT IS IMPORTANT TO PRESERVE THE ORDER OF vpot_site
    //IT MUST BE THE SAME AS USED TO SOLVE THE LINEAR SYSTEM OF EQUATIONS
    //YIELDING THE FACTORS TO CALCULATE a1, a2, ..., a7 
    for(int i = 0; i < Nbasis; i++){  
       pot  = &(FULL_ARRAY[count]); count += NH; 
       moms = &(FULL_ARRAY[count]); count += NH;
       spline_interpolation(ZH, pot, moms, NH, z1, &vpot, &dvd);
       vpot_site[i] = vpot; dvd_site[i] = dvd;
//       cout << i << " " << z1 << " " << vpot << "\n";
    }
    CALC_ATOMIC_UV(u0, v0, vpot_site, dvd_site, &vpot, &dvdu1, &dvdv1, &dvdz1);
  }
//  cout << z1 << " ";
  SET_TO_ORIGIN(&x1, &y1);
  double vref2; double dumx = dvdu1; double dumy = dvdv1; double dumz = dvdz1;
  REF_POT(x1, y1, z1, &vref, &dvdx1, &dvdy1, &dvdz11,0);
  vref2 = vref; dumx += dvdx1; dumy += dvdy1; dumz += dvdz11;
//  cout << vref << " "; 
  REF_POT(x1, y1, z1, &vref, &dvdx1, &dvdy1, &dvdz11,1);
  vref2 += vref; dumx += dvdx1; dumy += dvdy1; dumz += dvdz11;
// cout << vref << " "; 
  REF_POT(x1, y1, z1, &vref, &dvdx1, &dvdy1, &dvdz11,2);
  vref2 += vref;  dumx += dvdx1; dumy += dvdy1; dumz += dvdz11;
//  cout << vref << "\n";
   
  *vb1 = (vref2 + vpot)/27.211383;
//  *vb1 = vref2;
  dvd_cart[0] = -dumx/51.422065;
  dvd_cart[1] = -dumy/51.422065;
  dvd_cart[2] = -dumz/51.422065;
//   cout << pos_cart[0] << " " << pos_cart[1] << " " << pos_cart[2] << " " << vref2 + vpot << "\n";
};
void ATOMIC::CALC_ATOMIC_UV(double u0, double v0, double* vval, double* dvval, double* vpot, double* dvdu, double* dvdv, double* dvdz){
   //initialize BASIS FUNCTIONS AT THE POSITION u0 and v0
   //double u = gux*u0; double v = guy*v0;
   set_basis_val(u0, v0);
   double COEFF[Nbasis]; double sum = 0.; 
   double COEFF2[Nbasis];
   for(int i = 0; i < Nbasis; i++){ COEFF[i] = COEFF2[i] = 0.;}

   for(int i = 0; i < Nbasis; i++){
     for(int j = 0; j < Nbasis; j++){
       COEFF[i]  += vval[j]*INVERSE_Matrix[i*Nbasis+j]; 
       COEFF2[i] += dvval[j]*INVERSE_Matrix[i*Nbasis+j];
     }
   }

   for(int i = 0; i < Nbasis; i++) sum += COEFF[i]*get_basis(i);
   *vpot = sum;
   sum = 0.;
   for(int i = 0; i < Nbasis; i++) sum += COEFF[i]*get_dbasis_du(i);
   *dvdu = gux*sum;
   sum = 0.;
   for(int i = 0; i < Nbasis; i++) sum += COEFF[i]*get_dbasis_dv(i);
   *dvdv = guy*sum;
   sum = 0.;
   for(int i = 0; i < Nbasis; i++) sum += COEFF2[i]*get_basis(i);
   *dvdz = sum;
    
};