/*
    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. 30, 1997
 *
 * Calculates the gradient of the energy versus the electronic variables
 * Y (and also all the various routines needed to do this).
 *
 */

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

#include "header.h"
#include "parallel.h"

/*
 * Calculates the local part of the self-consistent potential.
 */
void calc_Vscloc(Elecvars &evars, Elecinfo &einfo, Symmetries &symm)
{
#ifdef DFT_PROFILING
  timerOn(25);   // Turn on calc_Vscloc timer
#endif // DFT_PROFILING
  int dir,dag; 

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

  switch (einfo.exc_opt) {
  case DFT_EXCORR_GGA_PBE:
    {
        RealSpaceScalarFieldColumn rtemp(n);
        CoeffSpaceScalarFieldColumn ctemp(d);

	// TAA: Test of gradient "D" functions (remove at cleanup)
//  	{
//  	  Basis *basis = n.basis;
//  	  Lattice *lattice = basis->basis_spec->lattice;

//  	  int Nx,Ny,Nz,Nx2,Ny2,Nz2,NyNz,NxNyNz;
//  	  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;
//  	  NxNyNz=basis->basis_spec->NxNyNz;
	  
//  	  RealSpaceScalarFieldColumn Dx(n);
	  
//  	  for (int i=-Nx2; i < Nx2; i++)
//  	    for (int j=-Ny2; j < Ny2; j++)
//  	      for (int k=-Nz2; k < Nz2; k++)
//  		{
//  		  int 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;
//  		  n.data.d[index] = cos(11*2*M_PI*j/((double) Nx))*lattice->R.m[0][0]/(2*M_PI);
//  		}
//  	  int dir,dag;
//  	  apply_D(dir=1,dag=0,n,Dx);
//  	  double tmp=0.;
//  	  FILE *fid=fopen("out","w");
//  	  for (int i=0; i<NxNyNz; i++)
//  	    {
//  	      tmp+=Dx.data.d[i].x*Dx.data.d[i].x+Dx.data.d[i].y*Dx.data.d[i].y;
//  	      fprintf(fid,"%20.10lf %20.10lf %20.10lf\n",n.data.d[i].x,Dx.data.d[i].x,Dx.data.d[i].y);
//  	    }
//  	  fclose(fid);
//  	  exit(1);
//  	}
	// End Test gradient function

        //Vscloc = Jdag(Vlocpot+O(d+J(exc_LDA(n,Dx))));
	RealSpaceScalarFieldColumn Dx(n),Dy(n),Dz(n);
  	apply_D(dir=0,dag=0,n,Dx);
  	apply_D(dir=1,dag=0,n,Dy);
  	apply_D(dir=2,dag=0,n,Dz);

        exc_GGA(n,Dx,Dy,Dz,rtemp);
        apply_J(rtemp,ctemp); /* TAA: Add in Hartree pot... */ ctemp+=d;
        apply_O(ctemp,ctemp); /* TAA: Add in loc pot... */ ctemp+=Vlocpot;
        apply_Jdag(ctemp, Vscloc);

        //Vscloc += pointwise_mult(Dnexc_GGA(n,Dx),Jdag(O(J(n))))
	//    + sum_i Ddag_i( pointwise_mult(DD_iexc_GGA(n,Dx),Jdag(O(J(n)))) )
        apply_J(n,ctemp); apply_O(ctemp,ctemp); 
        apply_Jdag(ctemp,rtemp);

        Vscloc+=pointwise_mult(Dnexc_GGA(n,Dx,Dy,Dz),rtemp);
	Vscloc+=Ddag(dir=0,pointwise_mult(DDexc_GGA(dir=0,n,Dx,Dy,Dz),rtemp));
	Vscloc+=Ddag(dir=1,pointwise_mult(DDexc_GGA(dir=1,n,Dx,Dy,Dz),rtemp));
	Vscloc+=Ddag(dir=2,pointwise_mult(DDexc_GGA(dir=2,n,Dx,Dy,Dz),rtemp));

        // and here is the place to call the symmetrizer Sdag to the potential
        symmetrize_n_dag(Vscloc, &symm);
        break;
    }
    case DFT_EXCORR_LDA:
      /*
       * Vscloc = Jdag(Vlocps) + Jdag(O(d)) + Jdag(O(J(exc(n)))) + 
       *            pointwise_mult(excprime(n),Jdag(O(J(n))));
       */
    {//ipd: why ptwise_mult- can't we overload *
        RealSpaceScalarFieldColumn rtemp(n);
        CoeffSpaceScalarFieldColumn ctemp(d);
        exc_LDA(n,rtemp);
        apply_J(rtemp,ctemp); ctemp += d;
        apply_O(ctemp,ctemp); ctemp += Vlocpot;	
        apply_Jdag(ctemp, Vscloc);
        //Vscloc = Jdag(Vlocpot+O(d+J(exc_LDA(n)))); 
        apply_J(n,ctemp); apply_O(ctemp,ctemp); 
        apply_Jdag(ctemp,rtemp);
        //Vscloc += pointwise_mult(excprime_LDA(n),Jdag(O(J(n))));
        Vscloc += pointwise_mult(excprime_LDA(n),rtemp);
        // and here is the place to call the symmetrizer Sdag to the potential
        symmetrize_n_dag(Vscloc, &symm);
        break; 
      }
      // Right now, we're only supporting PW LDA xc.  --mhe
      // And Teter's LSD.  --mhe
      
//    case DFT_EXCORR_GGA_PW91:
//      if (einfo.spintype == NOSPIN) {
//        Vscloc = Jdag(Vlocps+O(d+J(exc_LDA(n)+exprime_GC(n))));
//        Vscloc += pointwise_mult(excprime_LDA(n),Jdag(O(J(n))));
//      } else {
//        die("GGA_PW91 not implemented with spin\n");
//      }
//      break;

  case DFT_EXCORR_LSD_TETER:
    {
      RealSpaceScalarFieldColumn &Vscloc_up = evars.Vscloc_up;
      RealSpaceScalarFieldColumn &Vscloc_dn = evars.Vscloc_dn;

      RealSpaceScalarFieldColumn exc(Vscloc);
      RealSpaceScalarFieldColumn excprime_up(Vscloc);
      RealSpaceScalarFieldColumn excprime_dn(Vscloc);

      exc = exc_LSD_TETER(evars.n_up, evars.n_dn);
      excprime_LSD_TETER(evars.n_up, evars.n_dn, excprime_up, excprime_dn);

      Vscloc_up  = Jdag(Vlocpot);
      Vscloc_up += Jdag(O(d));
      Vscloc_up += Jdag(O(J(exc)));
      Vscloc_up += pointwise_mult(excprime_up,Jdag(O(J(n))));
      symmetrize_n_dag(Vscloc_up, &symm);
      
      Vscloc_dn  = Jdag(Vlocpot);
      Vscloc_dn += Jdag(O(d));
      Vscloc_dn += Jdag(O(J(exc)));
      Vscloc_dn += pointwise_mult(excprime_dn,Jdag(O(J(n))));
      symmetrize_n_dag(Vscloc_dn, &symm);
      break;
    }
    default:
      die("That exchange correlation is not supported.");
  }


#ifdef DFT_PROFILING
  timerOff(25);   // Turn off calc_Vscloc timer
#endif // DFT_PROFILING 
}

/*
 * Calls calc_Vscloc() [above] and then
 * calculates and diagonalizes the subspace Hamiltonian.
 */
void
calc_Hsub(Everything &everything)
{
  Elecvars &evars = everything.elecvars;
  Elecinfo &einfo = everything.elecinfo;
  Symmetries &symm = everything.symm;

  int q;
  /* Local workspaces */
  /*
   * to avoid allocate/deallocate many times, create for largest
   * col_length. 
   */
  //ipd: this won't work w/ the current structures
  // update it when we have scratch space manager
/*    int max_col_length = 0; */
/*    for (q=0; q < einfo.nstates; q++) */
/*      if (max_col_length < bs[q].C.col_length) */
/*        max_col_length = bs[q].C.col_length; */


  // if not doing band structure with reading the potential from file
  // we must compute it
  if( (everything.cntrl.band_structure_flag==FALSE)||(einfo.read_vscloc_flag==FALSE) )
    calc_Vscloc(evars,einfo,symm);
  
  
  for (q=0; q < einfo.nstates; q++) {
    BlochState &s = evars.states[q];
  
    ColumnBundle HspC(s.C);
    apply_Hsp(s.C, HspC, *s.Vscloc, everything);
    s.Hsub = (s.C)^HspC;
    s.Hsub.hermetian = 1;
    diagonalize_herm(s.Hsub_eigs,s.Hsub_evecs,s.Hsub,
                     s.Hsub.nr);
  }
}

/*
 * Call calc_Vscloc() [above] and then
 * calculates the energy gradient versus Y and B
 * and calculates and diagonalizes the subspace Hamiltonian.
 */
void
calc_elecgrad_and_Hsub(ColumnBundle **Ygrad, ComplexMatrix **Bgrad,
		       Everything &everything)
{

#ifdef DFT_PROFILING
  timerOn(10); // Turn on calc_elecgrad_... timer
#endif // DFT_PROFILING

  Elecinfo &einfo = everything.elecinfo;
  Elecvars &evars = everything.elecvars;
  Symmetries &symm = everything.symm;
  int q;

  /* Local workspaces */
  /*
   * to avoid allocate/deallocate many times, create for largest
   * col_length. 
   */
/*  int max_col_length = 0;
  for (q=0; q < einfo.nstates; q++)
    if (max_col_length < bs[q].C.col_length)
      max_col_length = bs[q].C.col_length;
*/
//  ColumnBundle tmp(bs[0].C.tot_ncols,max_col_length);
  calc_Vscloc(evars,einfo,symm);
  for (q=0; q < einfo.nstates; q++) {
    BlochState &s = evars.states[q];
    ColumnBundle tmp(s.C);
    /* Apply Hsp, calculate Hsub, and diagonalize Hsub */
    apply_Hsp(s.C,*Ygrad[q],*s.Vscloc, everything);
    s.Hsub = (s.C)^(*Ygrad[q]);
    s.Hsub.hermetian = 1;

    diagonalize_herm(s.Hsub_eigs,s.Hsub_evecs,s.Hsub,s.Hsub.nr);

    /* Calc Y gradient */
/* The code/mess below does: */
/*       HspC = Hsp(C[q],ioninfo) */
/*       Ygrad[q]  = Pbar(C[q],HspC*(F[q]*V[q]*Umhalf[q])); */
/*       Ygrad[q] += O(C[q]*V[q]*Q(Vdag[q]*(Hsub[q]*F[q]-F[q]*Hsub[q])*V[q],W[q],mu[q])); */
/*       Ygrad[q] *= w[q]; */
    if (einfo.subspace_rotation) {
        do_ColumnBundle_matrix_mult(*Ygrad[q],(s.F)*(s.V)*
                                    (s.Umhalf),tmp,0);
        apply_Pbar(s.C,tmp,*Ygrad[q]);
        do_ColumnBundle_matrix_mult(
            s.C,
            (s.V)*Q(herm_adjoint(s.V)*((s.Hsub)*(s.F)-
                                       (s.F)*(s.Hsub))*
                    (s.V),s.W,s.mu),
            tmp,0);
    } else {
      do_ColumnBundle_matrix_mult(*Ygrad[q],(s.F)*(s.Umhalf),tmp,0);
      apply_Pbar(s.C,tmp,*Ygrad[q]);
      do_ColumnBundle_matrix_mult(
        s.C,
        Q((s.Hsub)*(s.F)-(s.F)*(s.Hsub),s.W,
          s.mu),
        tmp,0);
    }
    // O(in,out) should work for this, right?
    //apply_O_inplace(tmp);
    O(tmp,tmp);
    *Ygrad[q] += tmp;
    *Ygrad[q] *= s.w;

    if (einfo.subspace_rotation) {
      /* Calc B gradient and symmetrize it */
      /* Bgrad[q] = R([Hsub[q],F[q]]) */
      *Bgrad[q] = R(((s.Hsub)*(s.F)-(s.F)*(s.Hsub)),s.Z,
                    s.beta);
      *Bgrad[q] += herm_adjoint(*Bgrad[q]);
      *Bgrad[q] *= (scalar)(0.5*(s.w));
    }
  }

#ifdef DFT_PROFILING
  timerOff(10); // Turn off calc_elecgrad_... timer
#endif // DFT_PROFILING
}

/*
 * Does what the above routine does plus preconditioning the gradient.
 * The preconditioning is diagonal inverse kinetic for the plane-wave basis.
 * Namely, the average kinetic energy per band is calculated 
 * (via ener->KE/einfo.nelectrons) and this (times 2.0) is used
 * as the roll-over value of the kinetic energy for the preconditioning.
 */
void
calc_elecgrad_pgrad_and_Hsub(ColumnBundle **Ygrad, ColumnBundle **pYgrad,
			     ComplexMatrix **Bgrad, ComplexMatrix **pBgrad,
			     Everything &everything)
{
#ifdef DFT_PROFILING
  timerOn(10); // Turn on calc_elecgrad_... timer
#endif // DFT_PROFILING

  Elecinfo &einfo = everything.elecinfo;
  Elecvars &evars = everything.elecvars;
  Symmetries &symm = everything.symm;

  int q;

  /* Local workspaces */
  /*
   * to avoid allocate/deallocate many times, create for largest
   * col_length. 
   */
  // can't do this anymore: problem? --csg
  // 
  //  int max_col_length = 0;
  //for (q=0; q < einfo.nstates; q++)
  //  if (max_col_length < s.C.col_length)
  //    max_col_length = s.C.col_length;
  //
  //  ColumnBundle tmp(bs[0].C.tot_ncols,max_col_length);

  /* Roll-over value for kinetic preconditioning:  2.0*average kinetic
   * energy per band. */
  // This is now computed inside K.
  //real KErollover;
  //KErollover = 2.0*ener.KE/einfo.nelectrons;

  calc_Vscloc(evars,einfo,symm);
  for (q=0; q < einfo.nstates; q++) {
    BlochState &s = evars.states[q];
    ColumnBundle tmp(s.C);

    /*
     * In this wonderful mess below, we calculate the subspace Hamiltonian,
     * the gradient, and the preconditioned gradient!  The formulae:
     *  C[q] = Y[q]*Umhalf[q]*Vdag[q]
     *  Hsub[q] = C[q]*Hsp(C[q])
     *  Ygrad[q] = w[q]*(g1 + g2)
     *  g1 = Pbar(Hsp(C[q]*F[q]*V[q]*Umhalf[q]))
     *     = Pbar(Hsp(Y[q]*Umhalf[q]*Vdag[q]*F[q]*V[q]*Umhalf[q]))
     *  g2 = O(C[q]*V[q]*Q(Vdag[q]*[Hsub[q],F[q]]*V[q]))
     *  pYgrad[q] = w[q]*(Pbar(precond(Pbar(Hsp(Y[q])))) + g2)
     */

    /* First we calculate Hsp(Y[q]), Hsp(Y[q])*Umhalf[q]*Vdag[q]=Hsp(C[q]), and Hsub[q] */

    // here pYgrad is the OUTPUT
    apply_Hsp(s.Y,*pYgrad[q],*s.Vscloc, everything);

    if (einfo.subspace_rotation)
      do_ColumnBundle_matrix_mult(*pYgrad[q],
                                  (s.Umhalf)*herm_adjoint(s.V), tmp,0);
    else
      do_ColumnBundle_matrix_mult(*pYgrad[q],(s.Umhalf),tmp,0);
      

    s.Hsub = (s.C)^tmp;

    /* Digonalize Hsub[q] */
    s.Hsub.hermetian = 1;
    diagonalize_herm(s.Hsub_eigs,s.Hsub_evecs,s.Hsub,
                     s.Hsub.nr);
    
    /* Now calculate g1 and the first term of pYgrad[q] */
    apply_Pbar(s.C,*pYgrad[q],tmp);
    if (einfo.subspace_rotation)
      do_ColumnBundle_matrix_mult(
        tmp,
        (s.Umhalf)*herm_adjoint(s.V)*(s.F)*(s.V)
        *(s.Umhalf),
        *Ygrad[q],0);
    else
      do_ColumnBundle_matrix_mult(tmp,(s.Umhalf)*(s.F)*
                                  (s.Umhalf),*Ygrad[q],0);
    //precond_inv_kinetic(tmp,KErollover);
    // ipd 09/14/02: Uncommenting the next 2 lines to the original state
    // looks like matt forgot about them :(
    K(tmp,everything);
    apply_Pbar(s.C,tmp,*pYgrad[q]);
    // Ygrad contains Pbar*Hsp*Y*Umhalf*F*Umhalf
    // tmp contains Pbar*Hsp*Y

    // Precondition only the first term in the gradient.
    // ipd 09/14/02: Commenting out the next 2 lines to the original state
    //*pYgrad[q] = *Ygrad[q];
    //K(*pYgrad[q],everything);

    /* Add in g2 = O*C*Q() term to Ygrad[q] and pYgrad[q] */
    if (einfo.subspace_rotation)
      do_ColumnBundle_matrix_mult(
        s.C,
        (s.V)*Q(herm_adjoint(s.V)*
                ((s.Hsub)*(s.F)-(s.F)*(s.Hsub))*
                (s.V),s.W,s.mu),
        tmp,0);
    else
      do_ColumnBundle_matrix_mult(s.C,
                                  Q((s.Hsub)*(s.F)-(s.F)*
                                    (s.Hsub),s.W,s.mu),
                                  tmp,0);
    //apply_O_inplace(tmp);
    O(tmp,tmp);
    *Ygrad[q] += tmp;
    // Use if we're only preconditioning the first term.
    *pYgrad[q] += tmp;
    // Precondition here.
    //*pYgrad[q] = *Ygrad[q]; 
    //K(*pYgrad[q],everything);

    /* scale by weights */
    *Ygrad[q] *= s.w;
    *pYgrad[q] *= s.w;
    
    /* Calc B gradient and symmetrize it */
    /* Bgrad[q] = R([Hsub[q],F[q]]) */
    if (einfo.subspace_rotation) {
      *Bgrad[q] = R(((s.Hsub)*(s.F)-(s.F)*(s.Hsub)),s.Z,
                    s.beta);
      *Bgrad[q] += herm_adjoint(*Bgrad[q]);
      *Bgrad[q] *= (0.5 * (s.w));
      *pBgrad[q] = ((scalar) einfo.pcond_for_Bgrad) * (*Bgrad[q]);
    }
  }
  
#ifdef DFT_PROFILING
  timerOff(10); // Turn off calc_elecgrad_... timer
#endif // DFT_PROFILING
  
}

/*
 * Calculates the energy gradient versus Y of the "energy" (objective)
 * function used in band structure calculations:
 *
 *         E_band = sum_q tr( C[q]^H*C[q] )
 */
void
calc_elecgrad_band_struct(ColumnBundle **Ygrad,
                          Everything &everything)
{
  
#ifdef DFT_PROFILING
  timerOn(10); // Turn on calc_elecgrad_... timer
#endif // DFT_PROFILING
  
  Elecinfo &einfo = everything.elecinfo;
  Elecvars &evars = everything.elecvars;
  
  int q;
  /* Local copies to make the formulae easier to read */
  //column_bundle *C = evars.C;
  //matrix *Umhalf = evars.Umhalf;
  
  /*
   * Local workspace:
   * to avoid allocate/deallocate many times, create for largest col_length. 
   */
  //int max_col_length = 0;
  //for (q=0; q < einfo.nstates; q++)
  //  if (max_col_length < s.C.col_length)
  //    max_col_length = s.C.col_length;
  //ColumnBundle tmp(bs[0].C.tot_ncols,max_col_length);
  
  for (q=0; q < einfo.nstates; q++){
    BlochState &s = evars.states[q];
  
    // manually adjust ColumnBundle length;
    //tmp.col_length = s.C.col_length;
    //copy_innards_ColumnBundle(&(s.C),&tmp);
    ColumnBundle tmp(s.C);

    /*
     * Calc Y gradient, with formula give by:
     *   Ygrad[q] = Pbar(H(C[q]))*U[q]^{-1/2}
     */
    apply_Hsp(s.C,*Ygrad[q],*s.Vscloc, everything);
    do_ColumnBundle_matrix_mult(*Ygrad[q],s.Umhalf,tmp,0);
    apply_Pbar(s.C,tmp,*Ygrad[q]);
    
  }
#ifdef DFT_PROFILING
  timerOff(10); // Turn off calc_elecgrad_... timer
#endif // DFT_PROFILING
}

/*
 * This routine calculates the gradient and preconditioned gradient
 * for band structure minimizations.  Preconditioning is done just
 * like above, except that since we don't have the kinetic energy
 * of the electrons explicitely, we *GUESS*:  the average kinetic
 * energy per electron for choosing the preconditioner roll-over point
 * is chosen to be 1 Hartree for now.
 */
void
calc_elecgrad_pgrad_band_struct(ColumnBundle **Ygrad,
                                ColumnBundle **pYgrad,
                                Everything &everything)
{
#ifdef DFT_PROFILING
  timerOn(10); // Turn on calc_elecgrad_... timer
#endif // DFT_PROFILING
  
  Elecinfo &einfo = everything.elecinfo;
  Elecvars &evars = everything.elecvars;

  int q;
  /* Local copies to make the formulae easier to read */
//    column_bundle *C = evars.C;
//    column_bundle *Y = evars.Y;
//    matrix *Umhalf = evars.Umhalf;
  
  /* Local workspaces:
   * to avoid allocate/deallocate many times create for largest col_length. */
  //int max_col_length = 0;
  //for (q=0; q < einfo.nstates; q++)
  //  if (max_col_length < bs[q].C.col_length)
  //   max_col_length = bs[q].C.col_length;
  //ColumnBundle tmp(bs[0].C.tot_ncols,max_col_length);

  /* Roll-over value for kinetic preconditioning:  2.0*average kinetic
   * energy per band. */
  // This is now done in K.
  //real kinetic_energy_per_electron = 1.0; // ener.KE/einfo.nelectrons
  //real KErollover;
  //KErollover = 2.0 * kinetic_energy_per_electron;
  for (q=0; q < einfo.nstates; q++) {
    BlochState &s = evars.states[q];
    // manually adjust ColumnBundle length;
    //tmp.col_length = bs[q].C.col_length;
    //copy_innards_ColumnBundle(&(bs[q].C),&tmp);
    ColumnBundle tmp(s.C);

    /* 
     * Compute gradient and preconditioned gradient.  Formulae:
     * Ygrad[q]  = Pbar(H(Y[q]))*U^-1
     * pYgrad[q] = Pbar(precondition(Pbar(H(Y[q]))))
     */
    apply_Hsp(s.Y,*pYgrad[q],*s.Vscloc, everything);
    apply_Pbar(s.C,*pYgrad[q],tmp);
    do_ColumnBundle_matrix_mult(tmp,(s.Umhalf)*(s.Umhalf),
                                *Ygrad[q],0);
    //precond_inv_kinetic(tmp,KErollover);
    K(tmp,everything);
    apply_Pbar(s.C,tmp,*pYgrad[q]);
  }
  
#ifdef DFT_PROFILING
  timerOff(10); // Turn off calc_elecgrad_... timer
#endif // DFT_PROFILING
  
}

/*
 * Prints the eigenenergies of the subspace Hamiltonian for each
 * state.
 */
void print_Hsub_eigs(Everything &everything,Output *out)
{
  Elecinfo &einfo = everything.elecinfo;
  Elecvars &evars = everything.elecvars;
  
  int q,i;
  
  out->printf("Band energies:\n");
  for (q=0; q < einfo.nstates; q++){
    out->printf("\nstate = %d   q_k = [ %lg %lg %lg ]   w = %lg",
                q,evars.states[q].qnum.kvec.v[0],
                evars.states[q].qnum.kvec.v[1],
                evars.states[q].qnum.kvec.v[2],
                evars.states[q].w);
    if (einfo.spintype == FREESPIN)
      out->printf("   spin = free\n");
    else
      out->printf("   spin = %d\n",evars.states[q].qnum.spin);
    
    out->printf("%4s  %13s  %13s  %13s\n",
		"band","filling   ","diag(Hsub) ","epsilon   ");
    out->printf("-------------------------------------------------\n");
    for (i=0; i < evars.states[q].Hsub.nr; i++)
      out->printf("%4d  %13.6le  %13.6le  %13.6le\n",
                  i,REAL(evars.states[q].F.c[i]),
                  REAL(evars.states[q].Hsub(i,i)),
                  evars.states[q].Hsub_eigs[i]);
  }
  out->flush();
}