/*********************************************************************************
* Copyright (C) 2013-2017 Alexey V. Akimov
*
* This file is distributed under the terms of the GNU General Public License
* as published by the Free Software Foundation, either version 2 of
* the License, or (at your option) any later version.
* See the file LICENSE in the root directory of this distribution
* or <http://www.gnu.org/licenses/>.
*
*********************************************************************************/
/**
  \file PW_methods1.cpp
  \brief The file implements methods for reading QE-generated wavefunctions in text format
*/

#include "PW.h"

/// liblibra namespace
namespace liblibra{

using namespace liblinalg;
using namespace libspecialfunctions;

/// libqobjects namespace
namespace libqobjects{


void wfc::aux_line2vec(string line,vector<double>& a){
// Speciall auxiliary function to conver line with 3 numbers to
// 3 distinct numbers in the array a
  vector<string> tmp; 
  split_line(line,tmp);
  if(a.size()>0) { a.clear(); }
  a = vector<double>(3,0.0);
  a[0] = atof(tmp[0].c_str());
  a[1] = atof(tmp[1].c_str());
  a[2] = atof(tmp[2].c_str());

}

void PW::QE_read_acsii_index(std::string filename){

  int verbose = 1;
  // Read all lines
  vector<std::string> A;
  int filesize = read_file(filename,1,A);
 
  // Get some information from the file
  int beg,end; beg = end = 0;
  int status = 0;
  size_t pos; string str;

  //-------------------- Set parameters to default values ------------
  nspin = -1;
  gamma_only = -1;
  natoms = -1; 
 
  //--------------------- Dimensions --------------------------------
  status = find_section(A,"<Dimensions>","</Dimensions>",0,filesize,beg,end);

  if(status){
    for(int i=beg;i<end;i++){
      pos = A[i].find("<Kpoints"); if(pos!=string::npos){ str = extract_s(A[i],"nktot"); nkpts = atoi(str.c_str()); }
      pos = A[i].find("<Kpoints"); if(pos!=string::npos){ str = extract_s(A[i],"nspin"); nspin = atoi(str.c_str()); }
      pos = A[i].find("<Wfc_grid");if(pos!=string::npos){ str = extract_s(A[i],"npwx"); npw = atoi(str.c_str()); }
      pos = A[i].find("<Bands"); if(pos!=string::npos){ str =  extract_s(A[i],"nbnd"); nbands = atoi(str.c_str()); }
      pos = A[i].find("<Gamma_tricks");if(pos!=string::npos){ str =  extract_s(A[i],"gamma_only"); if(str=="F"){ gamma_only=0; }else{ gamma_only=1;} }
      pos = A[i].find("<Atoms"); if(pos!=string::npos){ str =  extract_s(A[i],"natoms"); natoms = atoi(str.c_str()); }
    }// for i
  }// if status

  //------------------------ Cell ----------------------------------
  beg = end = 0;
  status = find_section(A,"<Cell","</Cell>",0,filesize,beg,end);

  if(status){
    for(int i=beg;i<end;i++){
      pos = A[i].find("<Cell"); if(pos!=string::npos){ str =  extract_s(A[i],"units"); cell_units = str; }
      pos = A[i].find("<Data"); if(pos!=string::npos){ str =  extract_s(A[i],"alat"); alat = atof(str.c_str()); }
      pos = A[i].find("<Data"); if(pos!=string::npos){ str =  extract_s(A[i],"omega"); omega = atof(str.c_str()); }
      pos = A[i].find("<Data"); if(pos!=string::npos){ str =  extract_s(A[i],"tpiba"); tpiba = atof(str.c_str()); }

      pos = A[i].find("<a1"); if(pos!=string::npos){ str =  extract_s(A[i],"xyz"); aux_line2vec(str,a1);  }
      pos = A[i].find("<a2"); if(pos!=string::npos){ str =  extract_s(A[i],"xyz"); aux_line2vec(str,a2); }
      pos = A[i].find("<a3"); if(pos!=string::npos){ str =  extract_s(A[i],"xyz"); aux_line2vec(str,a3); }
      pos = A[i].find("<b1"); if(pos!=string::npos){ str =  extract_s(A[i],"xyz"); aux_line2vec(str,b1); }
      pos = A[i].find("<b2"); if(pos!=string::npos){ str =  extract_s(A[i],"xyz"); aux_line2vec(str,b2); }
      pos = A[i].find("<b3"); if(pos!=string::npos){ str =  extract_s(A[i],"xyz"); aux_line2vec(str,b3); }

    }// for i
  }// if status

  //---------------------- Eigenvalues ----------------------------
  beg = end = 0;
  status = find_section(A,"<Eigenvalues","</Eigenvalues>",0,filesize,beg,end);

  if(status){
    for(int i=beg;i<end;i++){
      pos = A[i].find("<Eigenvalues"); if(pos!=string::npos){ str =  extract_s(A[i],"efermi"); efermi = atof(str.c_str()); }
      pos = A[i].find("<Eigenvalues"); if(pos!=string::npos){ energy_units =  extract_s(A[i],"units");  }
    }// for i

  }//status


  //======================= Almost done ============================
  // Before reading the eigenvalues we will allocate the memory
  // and set the corresponding flag so do not do this when reading wavefunctions
  // and grids

  if(nkpts>0 && nbands>0 && npw>0){

    if(is_allocated==0){ // Allocate all levels
      kpts = vector<K_point>(nkpts,K_point(nbands,npw));
      is_allocated = 3;
    }

  }// if >0 >0 >0


  // Now ready to read in the eigenvalues
  if(is_allocated==3){
    for(int k=0;k<nkpts;k++){
      int beg1,end1;
      // Use latest beg and end enclosing Eigenvalues section
      status = find_section(A,"<e."+int2str(k+1),"</e."+int2str(k+1),beg,end,beg1,end1);
 
      if(status){
        for(int i=0;i<nbands;i++){
          kpts[k].mo[i].energy = atof(A[beg1+1+i].c_str());
        }// for i
      }// if status
    }// for k
  }// is_allocated==3

}

