/*
    DFT++ is a density functional package developed by the research group
    of Professor Tomas Arias

    Copyright 1996-2003 Sohrab Ismail-Beigi

    This file is part of DFT++.

    DFT++ is free software; you can redistribute it and/or modify
    it 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.

    DFT++ is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with DFT++; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    Please see the file CREDITS for a list of authors.

    For academic users, we request that publications using results obtained with
    this software reference

    "New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi
    and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000).

    and, if using the wavelet basis, further reference

    "Multiresolution analysis of electronic structure: semicardinal and wavelet bases,"
    T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999).

    and 

    "Robust ab initio calculation of condensed matter: transparent convergence through
    semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and
    Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003).

    For your convenience, preprints of the above articles may be obtained from
    http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively.
*/

#include "header.h"

/////////// WL_BasisSpec //////////////////

// Default filename is "gridspec".
WL_BasisSpec::WL_BasisSpec()
{
  size_realgrid = 0;
  kpt_fold[0] = 0; kpt_fold[1] = 0; kpt_fold[2] = 0;
  basis_flag = 0;
}

// Setup function, called by the Everything class.
void WL_BasisSpec::setup(Everything &e)
{
  dft_log("\nSetting up wavelet BasisSpec.\n");

  // Set lattice to point to the global lattice class.
  lattice = &(e.lattice);
 
  FILE *fp;

  // Read grid structure from file 
  if ( (fp=fopen(e.basis_spec.filename,"r")) == NULL ) 
    die("Can't open gridspec file %s for input! \n",e.basis_spec.filename);

  // Modified readgridspec to not read a top-level scale.  This will be 
  // specified by the DFT++ Lattice class.
  readgridspec(fp,gridspec);
  fclose(fp);  

  // Now that we have the grid structure, get the scale ( L/#gridpts. )
//    gridspec[0].scale.x = lattice->latvec.m[0][0]/gridspec[0].dim.x;
//    gridspec[0].scale.y = lattice->latvec.m[1][1]/gridspec[0].dim.y;
//    gridspec[0].scale.z = lattice->latvec.m[2][2]/gridspec[0].dim.z;
  // This needs to be more general.  Lattice vectors are columns of latvec.
  vector3 a,b,c;
  a=lattice->latvec[0]; b=lattice->latvec[1]; c=lattice->latvec[2];
  gridspec[0].scale.x = abs(a)/gridspec[0].dim.x;
  gridspec[0].scale.y = abs(b)/gridspec[0].dim.y;
  gridspec[0].scale.z = abs(c)/gridspec[0].dim.z;  
  dft_log("\nTop level scale = [ %lf  %lf  %lf ] bohr\n",
          gridspec[0].scale.x, gridspec[0].scale.y, gridspec[0].scale.z);

  
  // compute the size of the grid
  int i=0;
  while ( gridspec[i].level >= 0 ) {
    size_realgrid+=(gridspec[i].dim.x*gridspec[i].dim.y*gridspec[i].dim.z);
    i++;
  }

  //now create the dense gridspec
  expandgridspec2x(gridspec,gridspec_dense);

  // if we are going to use dense grids - update the size
  if(use_dense)
    size_realgrid*=8;

  dft_log("Size of the realgrid %d\n", size_realgrid);
}

/////////// WL_Basis //////////////////////

WL_Basis::WL_Basis()
{
  basis_spec = NULL;
  qnum=NULL;

  nbasis = 0;
  lwork = 0;
}

void WL_Basis::setup_template(WL_BasisSpec &bs)
{
  // Set representation externally.
  // Default is realspace.
//  representation = REALSPACE;
  
  basis_spec = &bs;
  
  int i = 0; 
  int max_x = 0, max_y = 0, max_z = 0;
  nbasis = 0;
  
  // Calculate the total number of basis elements, as well as the size
  // of the largest grid.

  while(bs.gridspec[i].level >= 0) {
    nbasis += (bs.gridspec[i].dim.x)*(bs.gridspec[i].dim.y)*\
      (bs.gridspec[i].dim.z);

    if(bs.gridspec[i].dim.x > max_x)
      max_x = bs.gridspec[i].dim.x;
    if(bs.gridspec[i].dim.y > max_y)
      max_y = bs.gridspec[i].dim.y;
    if(bs.gridspec[i].dim.z > max_z)
      max_z = bs.gridspec[i].dim.z;
    i++;
  }

  dft_log("Maximum grid dimensions = [ %d  %d  %d ] \n",max_x,max_y,max_z);

  if(bs.use_dense)
    lwork = (2*max_x+3*HFL_O2)*(2*max_y+3*HFL_O2)*(2*max_z+3*HFL_O2);
  else
    lwork = (max_x+3*HFL_O2)*(max_y+3*HFL_O2)*(max_z+3*HFL_O2);
  
  return;
}

// ipd: Don't get confused here - the argument is just for namespace
// compatibility with the casee of a k-point dependent basis
// (planewaves)
void WL_Basis::setup(const vector3 &kvec)
{
  //qnum->kvec = kvec;

  return;
}