/****************************************************************************
 *
 * DFT++:  density functional package developed by
 *         the research group of Prof. Tomas Arias, MIT.
 *
 * Principal author: Sohrab Ismail-Beigi
 *
 * Modifications for MPI version: Kenneth P Esler,
 *                                Sohrab Ismail-Beigi, and
 *                                Tairan Wang.
 *
 * Modifications for LSD version: Jason A Cline,
 *                                Gabor Csanyi,
 *                                Tairan Wang, and
 *                                Evan Reed
 *
 * Modifications for lattice/Pulay forces: Gabor Csanyi and
 *                                         Sohrab Ismail-Beigi
 *
 * K-point symmetries:  Nicolaj Moll
 *
 * Ionic dynamics and relaxation:  Tairan Wang and
 *                                 Evan Reed
 * 
 * Copyright (C) 1996-2000 Sohrab Ismail-Beigi
 *
 * Permission granted to the academic community to 
 * freely distribute and use this software provided that this notice is
 * maintained intact, and all publications resulting from its use
 * cite the following reference:
 *
 * Sohrab Ismail-Beigi and T. A. Arias, "New Algebraic Formulation of
 * Density Functional Calculation", Computer Physics Communications,
 * volume 128, issue 1-2, pages 1-45, June 2000.
 *
 ****************************************************************************/

/*
 *     Gabor Csanyi           Jul 8, 2001
 *
 * Calculate electron density band by band, and dump it
 *
 */
#include <header.h>

void
dump_bands(Everything &e)
{

  /* Get local pointer to the relevant variables to make the formulae
   * legible */
  Elecinfo &einfo = e.elecinfo;
  Elecvars &evars = e.elecvars;
  Lattice &lattice = e.basis_spec.lattice;
  Symmetries &symm = e.symm;
  Ioninfo &ioninfo = e.ioninfo;


  diag_matrix *F = einfo.F;
  real *w = einfo.w;
  column_bundle *C = evars.C;
  RealSpaceColumn n(evars.n);   /* make a local copy of a density type thing */
  RealSpaceColumn n_up(evars.n_up);   /* template */
  RealSpaceColumn n_dn(evars.n_dn);   /* template */

  dft_log("\n>>>>>>>>>> Dumping the density band by band <<<<<<<<<<<<<<\n\n");
  dft_log_flush();

  dft_log("Calculating Hsub");  dft_log_flush();
  calc_UVCn(einfo,
	    evars,
	    symm,
	    lattice);
  calc_Hsub(e);
  dft_log("\n");  dft_log_flush();
  for (int q=0; q < einfo.nstates; q++)
    {
      do_column_bundle_matrix_mult(evars.C[q],
                                   evars.Hsub_evecs[q],
                                   evars.Y[q], 0);
      evars.C[q] = evars.Y[q];
    }
  
  int ncols = C[0].tot_ncols;
  dft_log("total col no=%d\n",ncols);
  dft_log_flush();
  char bandfilename[255];
  char bandfilenameup[255];
  char bandfilenamedn[255];
  
  if(einfo.spintype == NOSPIN)          // NOSPIN MODE
    { 
      for(int col = 0; col < ncols; col++){
        dft_log("computing col: %d ...",col);
        dft_log_flush();
        
        n.zero_out();
        for (int q=0; q < einfo.nstates; q++){
          n += ((scalar)w[q])*diagouterI_one_column(F[q],C[q],col);
        }
        dft_log("done\n",col);
        dft_log_flush();
        dft_log("writing to file...",col);
        dft_log_flush();
        sprintf(bandfilename, "banddens.%d", col);
        n.writea(bandfilename);
        dft_log("done\n",col);
        dft_log_flush();
        
      }
    }
  else if (einfo.spintype == ZSPIN)       // ZSPIN MODE
    {
      for(int col = 0; col < ncols; col++){
        dft_log("computing col: %d ...",col);
        dft_log_flush();
                
        n_up.zero_out();
        n_dn.zero_out();
        
	for (int q=0; q < einfo.nstates; q++){
	  if      (C[q].qnum->spin == 1)
	    n_up += ((scalar)w[q])*diagouterI_one_column(F[q],C[q],col);
	  else if (C[q].qnum->spin == -1)
	    n_dn += ((scalar)w[q])*diagouterI_one_column(F[q],C[q],col);
	  else
	    die("Bad spin designation in ZSPIN mode");
	}
        dft_log("done\n",col);
        dft_log_flush();
        dft_log("writing to file...",col);
        dft_log_flush();
            
	sprintf(bandfilenameup, "banddensup.%d", col);
	n_up.writea(bandfilenameup);
	sprintf(bandfilenamedn, "banddensdn.%d", col);
	n_dn.writea(bandfilenamedn);

        dft_log("done\n");
        dft_log_flush();
      }  
    }
  else 
    die("dump_bands is not implemented for this spintype\n");
}