/*
    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 all the various energy terms
 *
 */

/* $Id: calcener.cpp,v 1.11.2.26 2003/05/29 18:54:14 ivan Exp $ */

#include "header.h"

/*
 * Calclate kinetic energy:
 *
 *    KE = -0.5*sum_k { w[k]*trace(diag(C[k]^F[k]*L(C[k]))) }
 *
 */
void calc_KE(Everything &everything)
{
#ifdef DFT_PROFILING
  timerOn(18);   // Turn on calc_KE timer
#endif // DFT_PROFILING

  Elecinfo &einfo = everything.elecinfo;
  Elecvars &evars = everything.elecvars;
  Energies &ener = everything.energies;
    
  int q; // quantum state loop counter
  
  /* Create temporary workspace LC*/

  // now loop over states to accumulate the kinetic energy
  ener.KE = 0.0;
  for (q=0; q < einfo.nstates; q++){
    BlochState &s = evars.states[q];

    //ipd: for now create a copy - later use scratch space
    ColumnBundle LC(s.C);
    
    // ipd: ColumnBundle is much simpler than column_bundle
    // therefore we do not need to coppy the innards
    // copy_innards_column_bundle(s.C,&LC);
    
/* does:  ener.KE += REAL(sum_vector(diaginner(F[q],C[q],L(C[q]))))*w[q]; */
    L(s.C,LC);
    ener.KE +=  REAL(trace(diaginner(s.F,s.C,LC)))*s.w;
  }

  ener.KE *= -0.5;

#ifdef DFT_PROFILING
  timerOff(18);   // Turn off calc_KE timer
#endif // DFT_PROFILING
}

/*
 * Hartree energy.
 *
 * The routine assumes that the poisson equation has already been solved
 * and uses the electrostatic potential in evars.d.
 *
 */
void calc_EH(Everything &everything)
{
#ifdef DFT_PROFILING
  timerOn(21);   // Turn on calc_EH timer
#endif // DFT_PROFILING

  Elecvars &evars = everything.elecvars;

  RealSpaceScalarFieldColumn &n = evars.n;
  CoeffSpaceScalarFieldColumn &d = evars.d;

  //ipd: for WL the "^" operator changes n, just because
  //of our realspace conventions where we keep the ghosts
  //matched - try not to use it
  //everything.energies.EH = 0.5*REAL(n^Jdag(Obar(d)));
  
  everything.energies.EH = 0.5*REAL(J(n)^(Obar(d)));

#ifdef DFT_PROFILING
  timerOff(21);   // Turn off calc_EH timer
#endif // DFT_PROFILING
}

// Exchange correlation energy
void calc_Exc(Everything &everything)
{
#ifdef DFT_PROFILING
  timerOn(22);   // Turn on calc_Exc timer
#endif // DFT_PROFILING

  Elecinfo &einfo = everything.elecinfo;
  Elecvars &evars = everything.elecvars;
  Energies &ener = everything.energies;

  if (einfo.spintype == FREESPIN)
    die("FREESPIN is not implemented.");  

  RealSpaceScalarFieldColumn &n = evars.n;

  switch (einfo.exc_opt) {
  case DFT_EXCORR_GGA_PBE: {
    /* Evaluate gradient on real space grid */
    RealSpaceScalarFieldColumn Dx(n),Dy(n),Dz(n);
    int dir,dag;
    apply_D(dir=0,dag=0,n,Dx);
    apply_D(dir=1,dag=0,n,Dy);
    apply_D(dir=2,dag=0,n,Dz);

    ener.Exc = REAL(J(n)^(O(J(exc_GGA(n,Dx,Dy,Dz)))));
  }
  break;
  case DFT_EXCORR_LDA: {
    /* save on temps:
      exc(n,foo);
      J(foo, foo);
      O(foo,foo);
      Jdag(foo,foo);
    */
    ener.Exc = REAL(J(n)^(O(J(exc_LDA(n)))));
  }
  break;
  // For now, PW exchange-correlation is our only option. --mhe
    /*
  case DFT_EXCORR_GGA_PW91:
    if (einfo.spintype == NOSPIN) 
      ener.Exc = REAL(n^Jdag(O(J(exc_LDA(n)+ex_GC(n)))));
    else 
      die("GGA_PW91 not implemented with spin\n");
    break;
    */
  case DFT_EXCORR_LSD_TETER:
    {
      if (einfo.spintype == NOSPIN)
	die("Can't use LSD_TETER without spin\n");
      
      RealSpaceScalarFieldColumn &n_up = evars.n_up;
      RealSpaceScalarFieldColumn &n_dn = evars.n_dn;
      
      RealSpaceScalarFieldColumn exc(n);
      exc = exc_LSD_TETER(n_up, n_dn);
      ener.Exc = REAL(J(n)^(O(J(exc))));
      break;
    }
  default: 
    {
      die("This exchange correlation (%d) is not supported.",einfo.exc_opt);
      break;
    }
  }
  
  
#ifdef DFT_PROFILING
  timerOff(22);   // Turn off calc_Exc timer
#endif // DFT_PROFILING
}

// Ewald energy
void calc_Eewald(Everything &everything)
{
#ifdef DFT_PROFILING
  timerOn(24);   // Turn on calc_Eewald timer
#endif // DFT_PROFILING

  everything.energies.Eewald = Ewald(everything.ioninfo,
                                     everything.lattice);

#ifdef DFT_PROFILING
  timerOff(24);   // Turn on calc_Eewald timer
#endif // DFT_PROFILING
}

/*
 * Calculate local pseudopotential energy:
 *
 *    Eloc = real{  n^Jdag(Vlocps) } 
 *
 * THE ROUTINE USES THE LOCAL PSEUDOPOTENTIAL IN evars.Vlocps.
 * It does NOT recalculate evars.Vlocps.
 *
 */
void calc_Eloc(Everything &everything)
{
#ifdef DFT_PROFILING
  timerOn(19);   // Turn on calc_Eloc timer
#endif // DFT_PROFILING

  Elecvars &evars = everything.elecvars;

  RealSpaceScalarFieldColumn &n = evars.n;
  CoeffSpaceScalarFieldColumn &Vlocpot = evars.Vlocpot;

  everything.energies.Eloc = REAL(J(n)^(Vlocpot) );

#ifdef DFT_PROFILING
  timerOff(19);   // Turn off calc_Eloc timer
#endif // DFT_PROFILING
}


/* Total sum of all the above:  sets ener.Etot to the sum. */
void calc_Etot(Energies &ener)
{
  /* TAA ener.Etot  = ener.Exc; */
  ener.Etot  =
    ener.KE      +
    ener.Eloc    +
    ener.Enl     +
    ener.EH      + 
    ener.Exc     +
    ener.Ecore   + 
    ener.Eewald  +
    ener.Epulay;
}

/*
 * Calculates the Helmholtz free energy for the system :  F = Etot - T*S
 * Etot is assumed to be already calculated (ener.Etot).  The entropy S,
 * the temperature kT, the chemical potential mu, and the fillings and weights
 * are all inside of einfo.  This is used only in the case were we calculate 
 * fermi fillings (i.e. einfo.calc__fillings == 1).
 *
 * If w[q] are the state weights and f[q][i] is the i'th band's filling
 * at the q'th state, then (here fillings obey 0<f[q][i]<1 and
 * we assume we have a spin-compensated system, spintype = NOSPIN)
 *
 * T*S = -2.0*sum_q { w[q]*sum_i { f[q][i]*ln(f[q][i]) + 
 *                                 (1-f[q][i])*ln(1-f[q][i]) } }
 *
 */
void calc_F(Everything &e)
{
  real E,TS,fqi,kT;
  int q,i;

  if (e.elecinfo.calc_fillings == 0){
    // Don't die, just set F = Etot (==> S = 0, for constant fillings).
    // die("\n\ncalc_F() called with einfo.calcfillings == 0!!!\n\n");
    dft_log("For constant fillings, F = Etot.");
    e.energies.F = e.energies.Etot;
    return;
  }
    
  // Calculate T*S, and N 
  TS = (real)0.0;
  kT = e.elecinfo.kT;
  E = e.energies.Etot;
  for (q=0; q < e.elecinfo.nstates; q++)
    for (i=0; i < e.elecinfo.nbands; i++){
      fqi = REAL(e.elecvars.states[q].F.c[i])/2.0;
      if (fqi < 0.0 || fqi > 1.0)
	    die("\n\nf[%d][%d] is not in [0.0,1.0] in calc_G()!\n\n", q,i);

      // If fqi is not 0.0 or 1.0 (to one part in 10^30), then add its
      // contribution to TS 
      if ( fqi > 1.0e-30 && (1.0-fqi) > 1.0e-30)
        TS -= kT*2.0*e.elecvars.states[q].w*
          ( fqi*log(fqi) + (1.0-fqi)*log(1.0-fqi) );
    }
  e.energies.F = E - TS;
}


/*
 * Calculate the energy terms depenent on the electronic variables,
 * and recalculate the total energy.
 */
void calc_elec_dependent_energies(Everything &everything)
{
  Energies &ener = everything.energies;
  Elecinfo &einfo = everything.elecinfo;

  calc_KE(everything); 
  calc_Eloc(everything);
  calc_Enl(everything);
  calc_EH(everything);
  calc_Exc(everything);
  
  calc_Etot(ener); //add all the pieces

  if (einfo.calc_fillings == 1)
    calc_F(everything);
}

/*
 * Calculate the core and Ewald energies only and recalculate the total.
 */
//calc_ionic_energies(Everything &everything); //basis dependent
void calc_core_ewald_pulay_energies(Everything &everything)
{
  Energies &ener = everything.energies;
  Elecinfo &elecinfo = everything.elecinfo;

  calc_Ecore(everything); 
  calc_Eewald(everything);
  calc_Epulay(everything);
  calc_Etot(ener);
  
  if (elecinfo.calc_fillings == 1)
    calc_F(everything);
}

/*
 * Calculates ALL of the above energies
 */
//basis indep
void calc_all_energies(Everything &everything)
{
  calc_elec_dependent_energies(everything);
  calc_core_ewald_pulay_energies(everything);
}