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

// This is the generalized preconditioner, which wraps up 
// precond_inv_kinetic.
void K(ColumnBundle &c, Everything &e)
{
  /* Roll-over value for kinetic preconditioning:  2.0*average kinetic
   * energy per band. */
  real KErollover;
  if(e.energies.KE > 0)
    KErollover = 2.0*e.energies.KE/e.elecinfo.nelectrons;
  else
    KErollover = 2.0; // this takes care of the band structure case, where
                      // we don't know the KE per electron

  precond_inv_kinetic(c,KErollover);
}
  
/*
 * Multiply column_bundle by preconditioner.  In our case, we precond.
 * with the inverse kinetic part of Hsp.
 *
 * For each column of Y, we multiply the G'th component by
 * f(0.5*|k+G|^2/KErollover) (the kinetic energy) where
 * f(x) satisfies:
 *
 * f(x) = 1 + O(x^N)          for x << 1  (does nothing to low KE parts)
 *      = (1/x)*(1 + O(x^-N)) for x >> 1  (scale by 2*KErollover/|k+G|^2 for
 *                                               high KE components)
 *      = N/(N+1)             for x == 1
 *
 * f(x) = (1+x+x^2+x^3+...+x^(N-1))/(1+x+x^2+...+x^N)=(1-x^N)/(1-x^(N+1))
 * has this property, and it is used below.  Currently, N = 9 below.
 *
 * The function of KErollover is to determine the point where x=1, i.e.
 * the roll-over point for the function f(x).
 */
void precond_inv_kinetic(ColumnBundle &Y, real KErollover)
{
#ifdef DFT_PROFILING
  timerOn(33); // turn on precond timer
#endif // DFT_PROFILING

  int i,j;
  Basis *basis;
  Lattice *lattice;
  real kplusGx,kplusGy,kplusGz,kplusG2,kx,ky,kz;
  real GGTxx,GGTyy,GGTzz,GGTxy,GGTxz,GGTyz;
  real f,x;
  int nbasis,my_ncols;

  //  if (Y.basis == 0)
  //  die("K(ColumnBundle) called with Y.basis == 0\n");
  //basis = Y.basis;

  /* Let's go! */
  nbasis = Y.col_length;
  my_ncols = Y.my_ncols;
  basis = Y.col[0]->basis;
  lattice = Y.col[0]->basis_spec->lattice;
  kx = basis->qnum->kvec.v[0];
  ky = basis->qnum->kvec.v[1];
  kz = basis->qnum->kvec.v[2];
  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 (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    );
      x = 0.5*kplusG2/KErollover;
      f = 1.0+x*(1.0+x*(1.0+x*(1.0+x*(1.0+x*(1.0+x*(1.0+x*(1.0+x)))))));
      f = f/(1.0+x*f);

      // Hide ncols.
      for (i=0; i < my_ncols; i++)
        Y.col[i]->data.d[j] *= f;
    }

#ifdef DFT_PROFILING
  timerOff(33); // turn off precond timer
#endif // DFT_PROFILING
}