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

/*
 *     Sohrab Ismail-Beigi           Mar. 29, 1997
 *
 * Calculate the U matrices, the C column_bundles, and charge density n.
 *
 */

/* $Id: calcUVCn.cpp,v 1.13.2.12 2003/05/29 18:54:13 ivan Exp $ */

#include "header.h"

/*
 * Calculate the U matrices and their various powers:
 *
 *                       U[k] = Y[k]^O(Y[k])
 */
void calc_U(Elecinfo &einfo,Elecvars &evars)
{
#ifdef DFT_PROFILING
  timerOn(14);   // Turn on calc_U timer
#endif // DFT_PROFILING
 
  int q;
  
  for (q=0; q < einfo.nstates; q++) {
    BlochState &s=evars.states[q];
    ColumnBundle OY(s.Y);
    /* does:  U[k] = Y[k]^(O(Y[k])); */
    O(s.Y,OY);
    s.U = s.Y^OY;
    s.U.hermetian = 1;  // U is hermetian, so set the flag
    s.Umhalf = Uminusonehalf(s.U,s.W,s.mu);
  }
  
#ifdef DFT_PROFILING
  timerOff(14);   // Turn off calc_U timer
#endif // DFT_PROFILING
  
}

/*
 * Calculate the subspace rotation matrix  :  V[q] = exp(iB[q]) 
 */
void
calc_V(Elecinfo &einfo, Elecvars &evars)
{
  int q,i;

  /* If we aren't doing subspace rotations, then we should do nothing! */
  if (!einfo.subspace_rotation)
    return;
  
#ifdef DFT_PROFILING
  timerOn(41); // Calc V timer
#endif
  
  diag_matrix expibeta(evars.states[0].B.nr);//assuming everybody is the same
  // consider moving it inside the loop
  
  for (q=0; q < einfo.nstates; q++) {
    BlochState &s=evars.states[q];
    /* B[k] = Z[k]*beta[k]*Zdag[k] */
    diagonalize_herm(s.beta,s.Z,s.B,s.B.nr);
    /* V[k] = Z[k]*exp(i*beta[k])*Zdag[k] */
    for (i=0; i < s.B.nr; i++) {
#ifdef SCALAR_IS_COMPLEX
	  complex ibeta;
      
	  ibeta.x = 0.0;
	  ibeta.y = s.beta[i];
	  expibeta.c[i] = exp(ibeta);
#elif defined SCALAR_IS_REAL
#error "Don't know how to do subspace rotations for real"
#else
#error scalar is neither complex nor real
#endif
	}
    s.V = s.Z*(expibeta*herm_adjoint(s.Z));
  }

#ifdef DFT_PROFILING
  timerOff(41); // Calc V timer
#endif
}

/*
 * Calculate the orthonormal wavefunctions:
 *
 * state[q]: C = Y*Umhalf*Vdag
 */
void
calc_C(Elecinfo &einfo,Elecvars &evars)
{
#ifdef DFT_PROFILING
  timerOn(15);   // Turn on calc_C timer
#endif // DFT_PROFILING

  int q;
  //just loop over the bloch states and transform Y to C
  for (q=0; q < einfo.nstates; q++) {
    BlochState &s = evars.states[q];

    if (einfo.subspace_rotation)
      /*       C = Y*Umhalf; */
      do_ColumnBundle_matrix_mult(s.Y, s.Umhalf*herm_adjoint(s.V), s.C, 0);
    else
      do_ColumnBundle_matrix_mult(s.Y, s.Umhalf,s.C, 0);
  }
  
#ifdef DFT_PROFILING
  timerOff(15);   // Turn off calc_C timer
#endif // DFT_PROFILING
}

/* Combines U, V, and C */
void
calc_UVC(Elecinfo &einfo,Elecvars &evars)
{
  calc_U(einfo,evars);
  calc_V(einfo,evars);
  calc_C(einfo,evars);
}


/* Calculate electron density: n = sum_k { w[q]*diagouter(F[q],I*C[q]) } */
void
calc_n(Elecinfo &einfo, Elecvars &evars,
       Lattice &lattice, Symmetries &symm)
{
#ifdef DFT_PROFILING
  timerOn(16);   // Turn on calc_n timer
#endif // DFT_PROFILING


  /* Get local pointer to the relevant variables to make the formulae
   * legible */
  RealSpaceScalarFieldColumn &n = evars.n;   /* reference */
  RealSpaceScalarFieldColumn &n_up = evars.n_up;   /* reference */
  RealSpaceScalarFieldColumn &n_dn = evars.n_dn;   /* reference */
//  RealSpaceScalarFieldColumn &n_ud = evars.n_ud;   /* reference */
//  RealSpaceScalarFieldColumn &n_du = evars.n_du;   /* reference */

  //ipd: these guys are always present, no need of if sttatement
  n.zero_out();
  n_up.zero_out();
  n_dn.zero_out();


  //ipd: change the += below w/ multadds to save space
  for (int q=0; q < einfo.nstates; q++){
    BlochState &s = evars.states[q];
    switch(einfo.spintype){
      case NOSPIN:          // NOSPIN MODE
        n += ((scalar)s.w)*diagouterI(s.F,s.C); break;
      case ZSPIN:           // ZSPIN MODE
        switch (s.qnum.spin){
          case  1 : n_up += ((scalar)s.w)*diagouterI(s.F,s.C); break;
          case -1 : n_dn += ((scalar)s.w)*diagouterI(s.F,s.C); break;
          default : die("Bad spin designation in ZSPIN mode: %d", s.qnum.spin);
        }
        break;
      default:
        die("Invalid spintype: %d", einfo.spintype);
    }
  }
  switch(einfo.spintype){
    case NOSPIN:
      // symmetrize the total charge density.
      symmetrize_n(n, &symm);
      break;
    case ZSPIN:
        // symmetrize
      symmetrize_n(n_up, &symm);
      symmetrize_n(n_dn, &symm);
      n = n_up; n+= n_dn;
      break;
    default:
      die("Invalid spintype: %d", einfo.spintype);
  }
  
  // integral of n  should be  number of electrons per unit cell.
  /*ipd: it's more complicated for wavelets
    -- several choices:
    * do s^(JN) - may be expensive
    * add just the top level grids  - terrible
    * make a basis dependent function still have to deal w/ it's cost
    real sum = REAL(sum_vector(n));
    sum *= lattice.unit_cell_volume / n.n;
    if ( fabs(sum/einfo.nelectrons-1) > 1e-4 )
    {
    dft_log(DFT_ANAL_LOG,
    ">>>>Real space charge integrates to %g != %g (nelectrons)\n",
    sum, einfo.nelectrons);
    die("Quiting due to violation of charge conservation.");
    }
  */
  
#ifdef DFT_PROFILING
  timerOff(16);   // Turn off calc_n timer
#endif // DFT_PROFILING
}

// Combines U, V, C and n
void calc_UVCn(Elecinfo &einfo,Elecvars &evars,
	  Lattice &lattice,Symmetries &symm)
{
  calc_U(einfo,evars);
  calc_V(einfo,evars);
  calc_C(einfo,evars);
  calc_n(einfo,evars,lattice,symm);
}