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

void apply_O(const PW_ComplexColumn &c, PW_ComplexColumn &Oc)
{
#ifdef DFT_PROFILING
  timerOn(29); // turn on O timer
#endif // DFT_PROFILING

  if(&c!=&Oc)
    Oc = c;

#ifdef DFT_PROFILING
  timerOff(29); // turn off O timer
#endif // DFT_PROFILING
}

void apply_Obar(const PW_ComplexColumn &c, PW_ComplexColumn &Obarc)
{
   if(c.representation == REALSPACE)
    die("Trying to take Obar of realspace column!\n"); 

   if(&c != &Obarc)
     Obarc = c;

   Obarc.data.d[0].x = 0.0;

   return;
}

void apply_L(const PW_ComplexColumn &c, PW_ComplexColumn &Lc)
{
#ifdef DFT_PROFILING
  timerOn(31); // turn on L timer
#endif // DFT_PROFILING

  if(c.representation == REALSPACE)
    die("Trying to take L of realspace column!\n");
  Lc.representation=c.representation;
  
  int i,j,k, nbasis;
  int Nx,Ny,Nz,Nx2,Ny2,Nz2,NyNz, index;
  real kplusGx,kplusGy,kplusGz,kx,ky,kz;
  real GGTxx,GGTyy,GGTzz,GGTxy,GGTxz,GGTyz;
  real kplusG2, G2;

  Basis *basis = c.basis;
  Lattice *lattice = basis->basis_spec->lattice;

  if (basis == NULL)
    die("apply_L() called with c.basis == NULL\n");
  if ( c.length!=Lc.length )
    die("apply_L() called with different sizes for c and Lc\n");

  GGTxx = lattice->GGT.m[0][0];
  GGTyy = lattice->GGT.m[1][1];
  GGTzz = lattice->GGT.m[2][2];
  GGTxy = lattice->GGT.m[0][1];
  GGTxz = lattice->GGT.m[0][2];
  GGTyz = lattice->GGT.m[1][2];

  /* 
   * we need to decide if we have to compute G vectors, or if they 
   * are stored for us. 
   */

  if(basis->Gx == NULL || basis->Gy == NULL || basis->Gz == NULL){
    // we compute G vectors
    Nx = basis->basis_spec->Nx; Nx2 = Nx/2;
    Ny = basis->basis_spec->Ny; Ny2 = Ny/2;
    Nz = basis->basis_spec->Nz; Nz2 = Nz/2;
    NyNz = Ny*Nz;

    for (i=-Nx2; i < Nx2; i++)
      for (j=-Ny2; j < Ny2; j++)
	for (k=-Nz2; k < Nz2; k++)
	  {
	    index = 0;
	    if (k < 0) index += k+Nz;        else index += k;
	    if (j < 0) index += Nz*(j+Ny);   else index += Nz*j;
	    if (i < 0) index += NyNz*(i+Nx); else index += NyNz*i;
	    G2 = i*i*GGTxx +
	      j*j*GGTyy +
	      k*k*GGTzz +
	      2*(i*j*GGTxy +
		 i*k*GGTxz +
		 j*k*GGTyz);
	    Lc.data.d[index] = -G2*c.data.d[index];
	  }
  } else {
    // G vectors are stored

    // make sure we have a quantum number
    if(basis->qnum == NULL)
      die("In apply_L() G vectors are not stored and qnum == NULL\n");

    kx = basis->qnum->kvec.v[0];
    ky = basis->qnum->kvec.v[1];
    kz = basis->qnum->kvec.v[2];
    nbasis = basis->nbasis;

    for (j=0; j < nbasis; j++){
      kplusGx = kx + basis->Gx[j];
      kplusGy = ky + basis->Gy[j];
      kplusGz = kz + basis->Gz[j];
      kplusG2 = (kplusGx*kplusGx*GGTxx +
		 kplusGy*kplusGy*GGTyy +
		 kplusGz*kplusGz*GGTzz +
		 2*(kplusGx*kplusGy*GGTxy +
		    kplusGx*kplusGz*GGTxz +
		    kplusGy*kplusGz*GGTyz    ));
      
        Lc.data.d[j] = -kplusG2*c.data.d[j];
    }
  }

#ifdef DFT_PROFILING
  timerOff(31); // turn off L timer
#endif // DFT_PROFILING

}

/*
 * The inverse of the Laplacian on vectors:  the is obviously modulo
 * G=0 where there are infinities.  For planewaves, we have -1/|G|^2.
 */
void apply_invL(const PW_ComplexColumn &c, PW_ComplexColumn &invLc)
{
#ifdef DFT_PROFILING
  timerOn(32); // turn on invL timer
#endif // DFT_PROFILING

  int Nx,Ny,Nz,Nx2,Ny2,Nz2,NyNz;
  int i,j,k,index;
  real G2,GGTxx,GGTxy,GGTxz,GGTyy,GGTyz,GGTzz;

  BasisSpec *spec;
  Lattice *lattice;

  if (c.basis == 0)
    die("invL(vector) called with v.basis==0\n");
  if (c.length != c.basis_spec->NxNyNz)
    die("invL(vector) called with vector.n != FFT box size\n");
  if ((c.basis_spec->Nx%2)!=0 || (c.basis_spec->Ny%2)!=0 || 
      (c.basis_spec->Nz%2)!=0)
    die("invL(vector) called with basis->Nx,Ny,or Nz not a multiple of 2\n");

  spec=c.basis_spec;
  lattice = spec->lattice;
  Nx = spec->Nx; Nx2 = Nx/2;
  Ny = spec->Ny; Ny2 = Ny/2;
  Nz = spec->Nz; Nz2 = Nz/2;
  NyNz = Ny*Nz;
  GGTxx = lattice->GGT.m[0][0];
  GGTyy = lattice->GGT.m[1][1];
  GGTzz = lattice->GGT.m[2][2];
  GGTxy = lattice->GGT.m[0][1];
  GGTxz = lattice->GGT.m[0][2];
  GGTyz = lattice->GGT.m[1][2];
  for (i=-Nx2; i < Nx2; i++)
    for (j=-Ny2; j < Ny2; j++)
      for (k=-Nz2; k < Nz2; k++)
	{
	  /* G = 0 case */
	  if (i==0 && j==0 && k==0)
	    invLc.data.d[0] = 0.0;
	  else
	    {
	      index = 0;
	      if (k < 0) index += k+Nz;        else index += k;
	      if (j < 0) index += Nz*(j+Ny);   else index += Nz*j;
	      if (i < 0) index += NyNz*(i+Nx); else index += NyNz*i;
	      G2 = i*i*GGTxx +
		   j*j*GGTyy +
		   k*k*GGTzz +
		   2*(i*j*GGTxy +
		      i*k*GGTxz +
		      j*k*GGTyz);
	      invLc.data.d[index] = -c.data.d[index]/G2;
	    }
	}

#ifdef DFT_PROFILING
  timerOff(32); // turn off invL timer
#endif // DFT_PROFILING
  return;
}