/*
    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.
*/

/*************************************************************************
 *
 * Defines the PlaneWave operations (mostly for the Solid State case)
 *
 *************************************************************************/

#include "header.h"

//constructors
PW_ComplexColumn::PW_ComplexColumn()
{
  fat = TRUE;
  length = 0;
  basis_spec = NULL;
  basis = NULL;
  representation = REALSPACE; //default

  embedding = SPHERE;
}

PW_ComplexColumn::PW_ComplexColumn(PW_Basis *b, char s[10])
{
  if(!strcmp(s,"fat"))
    fat = TRUE;
  else if (!strcmp(s,"slim"))
    fat = FALSE;
  else
    die("PW_Column(): Identity crisis - don't know if i'm slim or fat\n");

  basis = b;

  representation = COEFFSPACE; //default

  if( b != NULL)
    basis_spec = b->basis_spec;
}

PW_ComplexColumn::PW_ComplexColumn(PW_Basis *b, int space, char s[10])
{
  if(!strcmp(s,"fat"))
    fat = TRUE;
  else if (!strcmp(s,"slim"))
    fat = FALSE;
  else
    die("PW_Column(): Identity crisis - don't know if i'm slim or fat\n");

  basis = b;
  representation = space;

  if( b != NULL)
    basis_spec = b->basis_spec;
}

//copy constructor w/ extra argument
PW_ComplexColumn::PW_ComplexColumn(const PW_ComplexColumn &col,
                                   char s[10])
{
  if(!strcmp(s,"fat"))
    fat = TRUE;
  else if (!strcmp(s,"slim"))
    fat = FALSE;
  else
    die("PW_Column(): Identity crisis - don't know if i'm slim or fat\n");

  length = col.length;
  basis_spec = col.basis_spec;
  basis = col.basis;
  representation = col.representation; 
  embedding = col.embedding;

  if(fat) {
    data.init(length);
    data = col.data;
  }
  else
    die("You can't copy a slim column, fatso!\n");
    
}

/*
//copy constructor
PW_ComplexColumn::PW_ComplexColumn(const PW_ComplexColumn &col)
{
  fat = col.fat;
  length = col.length;
  basis_spec = col.basis_spec;
  basis = col.basis;
  representation = col.representation; 
  embedding = col.embedding;

  if(fat) {
    data.init(length);
    data = col.data;
  }
  else
    die("You can't copy a slim column, fatso!\n");
}
*/
// destructor
PW_ComplexColumn::~PW_ComplexColumn()
{
  if(fat)
    data.free();
}

void
PW_ComplexColumn::init_wavefunction(PW_Basis &b)
{
  length = b.nbasis;
  basis_spec = b.basis_spec;
  basis = &b;
  representation = COEFFSPACE; 
  embedding = SPHERE;
}


void
PW_ComplexColumn::init_scalarfield(PW_Basis &b, int space)
{
  basis = &b;
  basis_spec = b.basis_spec;
  length = basis_spec->NxNyNz;
  embedding = BOX;

  representation = space;

  if(fat)
    data.init(length);
  else
    die("You cannot initialize a slim column to contain a scalar field\n");
}



void 
PW_ComplexColumn::print()
{
  int j;

  for (j = 0; j < length; j++)
    dft_log("%f\t%f\n", data.d[j].x, data.d[j].y);

}

char *
PW_ComplexColumn::getrepresentation()
{
  static char realspacename[] = "RealSpace";
  static char coeffspacename[] = "CoeffSpace";

  switch(representation){
    case REALSPACE:
      return realspacename;
    case COEFFSPACE:
      return coeffspacename;
    default:
      die("Unknown representation: %d\n", representation);
  }

  return NULL;
}


char *
PW_ComplexColumn::getembedding()
{
  static char spherename[] = "Sphere";
  static char boxname[] = "Box";

  switch(embedding){
    case SPHERE:
      return spherename;
    case BOX:
      return boxname;
  default:
    die("Unknown embedding: %d\n", embedding);
  }

  return NULL;
}


// LINEAR ALGEBRA  - just call the array algebra
void
PW_ComplexColumn::operator=(const PW_ComplexColumn &col)
{
  if(length != col.length)
    die("length != col.length in PW_ComplexColumn::operator=(): %d %d\n", length, col.length);
  if(embedding != col.embedding)
    die("embedding != col.embedding in PW_ComplexColumn::operator=()\n");

  representation = col.representation;

  data = col.data;

}


void
PW_ComplexColumn::operator+=(const PW_ComplexColumn &col)
{
  if(length != col.length)
    die("length != col.length in PW_ComplexColumn::operator+=(): %d %d\n", length, col.length);
  if(representation != col.representation)
    die("representation != col.repr. in PW_ComplexColumn::operator+=()\n");
  if(embedding != col.embedding)
    die("embedding != col.embedding in PW_ComplexColumn::operator+=()\n");


  data += col.data;
}

void
PW_ComplexColumn::operator*=(const PW_ComplexColumn &col)
{
  if(length != col.length)
      die("length != col.length in PW_ComplexColumn::operator*=() %d %d\n", length, col.length);
  if(representation != col.representation)
    die("representation != col.repr. in PW_ComplexColumn::operator*=()\n");
  if(embedding != col.embedding)
    die("embedding != col.embedding in PW_ComplexColumn::operator*=()\n");


  data *= col.data;
}


void
PW_ComplexColumn::operator-=(const PW_ComplexColumn &col)
{
  if(length != col.length)
    die("length != col.length in PW_ComplexColumn::operator-=()%d %d\n", length, col.length);
  if(representation != col.representation)
    die("representation != col.repr. in PW_ComplexColumn::operator-=()\n");
  if(embedding != col.embedding)
    die("embedding != col.embedding in PW_ComplexColumn::operator-=()\n");


  data -= col.data;
}

void
PW_ComplexColumn::operator*=(const real r)
{
  data*=r;
}

void
PW_ComplexColumn::operator*=(const complex c)
{
  data*=c;
}


complex
PW_ComplexColumn::operator^(const PW_ComplexColumn &col)
{
  return dot(this->data, col.data);
}



// BASIS SPECIFIC


/* Fill the column with random gaussian distributed numbers
 * with zero mean and std. dev. = 1/(1+((0.5*G^2)/0.75)^6), where 0.5*G^2
 * is the kinetic energy of the basis function in question.
 * The reason for this strange formula is that we want high-energy
 * plane-waves to have lower initial weight than low-energy one.  The
 * cutoff of 0.75 is because following the Jones-scheme business for
 * the 2-atom Si cell, we want to keep the (0 0 0), (1 1 1), and (2 0 0)
 * vectors (units are 2*pi/a).  The energy of the (2 0 0) in Hartrees
 * and for a=5.43 Angstorms is 0.7499 (miracle!?!?), and the next
 * highest planewave is (2 2 0) with energy 1.4998.  So I just chose
 * some function which turns over at 0.75 relatively sharply.
 *
 * -- SIB
 */

void
PW_ComplexColumn::randomize()
{
  int j;
  real std,KE,t, ttt;
  vector3 kplusG;

  // ipd: hmm, sth's fishy here -should have it work for anything
  // this basically limits us to work on spheres
  //if (basis->nbasis != length)
  //  die("\nPW_ComplexColumn::randomize() nbasis != length!!!\n\n");

  switch (representation){
    case REALSPACE:
      //just put random stuff in there
      data.randomize();
      break;
    case COEFFSPACE:
      if(embedding==SPHERE)
        for(j=0; j < basis->nbasis; j++){
          kplusG.v[0] = (real)basis->Gx[j] + basis->qnum->kvec.v[0];
          kplusG.v[1] = (real)basis->Gy[j] + basis->qnum->kvec.v[1];
          kplusG.v[2] = (real)basis->Gz[j] + basis->qnum->kvec.v[2];
          KE = 0.5*(kplusG*(basis->basis_spec->lattice->GGT*kplusG));
          t = KE/0.75;
          //std = 1.0/(1.0+t*t*t*t*t*t);
          ttt = t*t*t;
          t = ttt*ttt;
          std = 1.0/(1.0+t);
          
          data.d[j].x = gauss(std);
          data.d[j].y = gauss(std);
        }
      else //embedding == BOX
        //i'll put uniform noise for now
        data.randomize();
        //gabor, looks we need this just for the tests - you can put
        //k-dep if you want, i'm leaving to you the resurection ...
/*          die("Don't know how to randomize %s %s\n", */
/*              getrepresentation(), getembedding()); */
      break;
    default:
      die("unknown representation %s in randomize\n", getrepresentation());
  }
  
}


void
PW_ComplexColumn::setmodones()
{
  int ii;
  real theta = rand(); // just a random phase
  complex one(cos(theta), sin(theta));

  for(ii=0; ii<length; ii++)
    data.d[ii] = one;

  //PW_ComplexColumn out(*this); // create a copy column
  if(this->representation == COEFFSPACE)
    apply_J(*this, *this);

  return;
}







void
PW_ComplexColumn::IdentityMap(PW_ComplexColumn &col) const
{
  if(this != &col)
    col = *this;
}



void 
PW_ComplexColumn::SphereToBoxMap(PW_ComplexColumn &col) const
{
  int *index = this->basis->index;
  int j;

  //ipd: considert col.data.zero_out();
  // Fill in part of G-space taken up by the basis with correct values 
  for (j=0; j < col.basis->basis_spec->NxNyNz; j++){
//    dft_log("zeroing %d\n", j);
    col.data.d[j] = 0.0;
  }

  for (j=0; j < basis->nbasis; j++){
    //dft_log("mapping %d\n", j);
    col.data.d[index[j]] = this->data.d[j];
  }

}

void 
PW_ComplexColumn::BoxToSphereMap(PW_ComplexColumn &col) const 
{
  int *index = col.basis->index;
  int j;

  // Fill in part of G-space taken up by the basis with correct values 
  for (j=0; j < col.basis->nbasis; j++){
    //dft_log("UNmapping %d\n", j);
    col.data.d[j] = this->data.d[index[j]];
  }

}

//ipd: adding the boxtospheremap and the mapflag functionality
void
PW_ComplexColumn::map(PW_ComplexColumn &mappedCol) const
{
  switch(embedding){
  case BOX:
    switch(mappedCol.embedding){
    case BOX:
      this->IdentityMap(mappedCol);
      break;
    case SPHERE:
      this->BoxToSphereMap(mappedCol);
      break;
    default:
      die("Invalid embedding in PW_ComplexColumn::map() : %d\n", mappedCol.embedding);
    }
    break;
  case SPHERE:
    switch(mappedCol.embedding){
    case SPHERE:
      this->IdentityMap(mappedCol);
      break;
    case BOX:
      this->SphereToBoxMap(mappedCol);
      break;
    default:
      die("Invalid embedding in PW_ComplexColumn::map() : %d\n", mappedCol.embedding);
    }
    break;
  default:
    die("Invalid map in PW_ComplexColumn::map() : %d\n", embedding);
  }
}

void
PW_ComplexColumn::read(char *filename)
{
  data.read(filename);
}


void
PW_ComplexColumn::write(char *filename)
{
  data.write(filename);
}




// adds the square of the absolute value
void
add_scale_abs2(const complex &c, const  PW_ComplexColumn & in,
                    PW_ComplexColumn &out)
{
  add_scale_abs2(c, in.data, out.data);
}



//returns the absolute squares of the input column
PW_ComplexColumn abs2(PW_ComplexColumn &in)
{
  int i;
  PW_ComplexColumn out(in);

  for(i=0; i< in.length; i++)
    out.data.d[i]=abs2(in.data.d[i]);

  return out;  
}

PW_ComplexColumn
pointwise_mult(const PW_ComplexColumn &a,
	       const PW_ComplexColumn &b)
{
   if(a.length != b.length)
    die("a.length != b.length in pointwise_mult!\n");

  PW_ComplexColumn ab(a);
  int i;

  for(i=0; i<a.length; i++)
    ab.data.d[i] *= b.data.d[i];

  return ab; 
}

// UNUSED CODE



/* Fill the column_bundle with random gaussian distributed numbers
 * with zero mean and std. dev. = 1/(1+((0.5*G^2)/0.75)^6), where 0.5*G^2
 * is the kinetic energy of the basis function in question.
 * The reason for this strange formula is that we want high-energy
 * plane-waves to have lower initial weight than low-energy one.  The
 * cutoff of 0.75 is because following the Jones-scheme business for
 * the 2-atom Si cell, we want to keep the (0 0 0), (1 1 1), and (2 0 0)
 * vectors (units are 2*pi/a).  The energy of the (2 0 0) in Hartrees
 * and for a=5.43 Angstorms is 0.7499 (miracle!?!?), and the next
 * highest planewave is (2 2 0) with energy 1.4998.  So I just chose
 * some function which turns over at 0.75 relatively sharply.
 */


/*  void */
/*  PW_Column<RealArray>::randomize(void) */
/*  { */
/*    register int i,j; */
/*    register real std,KE,t,ttt; */
/*    vector3 kplusG; */

/*    if (basis->nbasis != length) */
/*      die("\nPW_Column<RA>::randomize() nbasis != length!!!\n\n"); */

/*    for(j=0; j < basis->nbasis; j++){ */
/*        kplusG.v[0] = (real)basis->Gx[j] + qnum->kvec.v[0]; */
/*        kplusG.v[1] = (real)basis->Gy[j] + qnum->kvec.v[1]; */
/*        kplusG.v[2] = (real)basis->Gz[j] + qnum->kvec.v[2]; */
/*        KE = 0.5*(kplusG*(basis->lattice->GGT*kplusG)); */
/*        t = KE/0.75; */
/*        //std = 1.0/(1.0+t*t*t*t*t*t); */
/*        ttt = t*t*t; */
/*        t = ttt*ttt; */
/*        std = 1.0/(1.0+t); */

/*        this->data.d[j] = gauss(std); */
/*    } */
/*  } */





// this will be different, but let's just leave it like this for now,
// until we figure out the gamma point case
/*  void */
/*  PW_Column<RealArray>::init_wavefunction(PW_Basis &b) */
/*  { */
/*    basis = &b; */
/*    basis_spec = b.basis_spec; */
/*    map = &(PW_Column<RealArray>::SphereToBoxMap); */
/*  } */