/*
    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     August 2 1999
//
// This is the main DFT++ executable.
// It handles electronic minimizations, calculation of forces, and
// ionic dynamics.  Ionic relaxations are not working yet.
// It has been cleaned up a bit and the structure should be clear.
// The basic idea is to do any MPI setup if needed, then
// to try to read the input file, and then to do any minimizations/force/
// dynamics as needed, and finally to dump the current state and exit.
//

/* $Id: dft-two.cpp,v 1.1.2.92 2003/05/29 18:54:16 ivan Exp $ */

#include "header.h"
#include "parallel.h"
#include <time.h>

// Default input filename
#define DEFAULT_INPUTFILE    "dft.in"

// This report level reports all the basic messages we might
// be interested to the log file.  More or less anal options
// can be found in header.h
#define REPORT_LEVEL       DFT_BASIC_LOG

// Prints a banner
void print_banner(char **argv,char *input_filename);

// Parses the command line arguments to main().  See down below.
void parse_command_line(int argc, char **argv, 
			char *input_filename,
			char *log_filename);

void optest(Everything &e);

void Etest(Everything &e);

// HACK ALERT!!!
// This variable stores the square magnitude of Ygrad globally, so that the
// wavelet invL routing can access it.
double global_abs2Yg;

/***************************************************************************
 *                                                                         *
 * Main Program:                                                           *
 *                                                                         *
 ***************************************************************************/
int
main(int argc,char**argv)
{
  // Do any system dependent setup if needed
  System::GlobalInit(&argc,&argv);  

#ifdef DFT_PROFILING
  // Activate timers and counters
  timerActivateAll();
  counterActivateAll();
  // Turn initialization timer on
  timerOn(0); // Initialization timer on
#endif

  // The variables needed to do the work.  Contains inside everything
  // we will need!
  Everything everything;

  /**********************************************************
   **                                                      **
   ** Reading command line, setting up logfile.            **
   **                                                      **
   **********************************************************/
  // initialize input/output files to defaults.
  char input_filename[DFT_FILENAME_LEN] = DEFAULT_INPUTFILE;
  char log_filename[DFT_FILENAME_LEN] = "";

  // Parse command line and get input and log filenames.
  parse_command_line(argc,argv,input_filename,log_filename);

  // Open the logfile and make the system logger point to it
  Output logfile(REPORT_LEVEL,log_filename);
  System::global_log = &logfile;

  /**********************************************************
   ** Setup section: read and parser the input file, and   **
   ** the preform the various setup routines to prepare    **
   ** everything.                                          **
   **********************************************************/
  // Try to open the input file
  dft_text_FILE *inputfile = dft_text_fopen(input_filename,"r");
  if (inputfile == NULL)
    die("\n%s:  can't read '%s'.  Aborting.\n\n",argv[0], input_filename);

  // Print a welcome banner with useful information
  print_banner(argv,input_filename);

  // expand all includes in the input file.
  inputfile->expand_include();

  // parse input file and do the setup
  parse(inputfile,everything);
  dft_text_fclose(inputfile);

  // ipd: hack to check it works - move it in everything.setup
  // and put the actual variable names in there
//   System::setup_scratch(3, 50, sizeof(complex));
//   System::setup_scratch(3, 25, sizeof(complex));
//   System::setup_scratch(3, 30, sizeof(complex));

  // initial setup
  everything.setup();

#ifdef DFT_PROFILING
  timerOff(0);  // turn off initialization timer.
  timerOn(1);   // turn on total computation timer.
#endif

  //optest(everything); exit(1);

  // Calculate the local pseudopotential vector if not reading it from file
  //
  // NOTE: this routine is basis dependent, and so
  // has to be supplied by the basis developers
  //
  // if doing band structure we supply the self consistent potential
  // which would include Vloc from the ions.
  // change this once starting from charge densities is included.
  if(everything.cntrl.band_structure_flag != TRUE)
    Vlocpot(everything);  

  // calculate electron independent energies just once
  calc_core_ewald_pulay_energies(everything);

            
  
  //TEST the symemtries
//  symmetrize_n(everything.elecvars.Vlocpot, &everything.symm);
//  everything.elecvars.Vlocpot.write("Vs");


  // here we may try to 
  if(everything.cntrl.finite_diff_test_flag) {
    finitedifftest(everything);
    die("FINITE DIFF TEST DONE, EXITTING\n");
  }

    
  //calc_all_energies(everything);
  minimize_elec(everything);

  //now do the forces - no iondynamics loop yet
  if ( everything.cntrl.calculate_forces_flag ) {
    calc_UVCn_d(everything);
    // Calculate and report the forces
    calc_ionic_forces(everything);
    everything.ioninfo.print_force(dft_global_log, everything);
  }
  
  // Setup the signal handlers: turned off for now...
  // the signal setup is disable right now. only a skeleton is left.
  // setup_signals(&elecinfo,&elecvars);


#ifdef DFT_PROFILING
  timerOff(1); // turn off total computation timer.
#endif

  
  // Write out final electronic state to disk
  dft_log("\nDone!  Dumping final variables:\n\n");
  dft_log_flush();
  dump_and_stamp_variables(everything);


#ifdef DFT_PROFILING
  // Print a final timer and counter report (with some very large
  // iteration number, which should be our last!
  timer_counter_Report(9999);
#endif

  // End any system dependent stuff
  System::GlobalFinalize();

  return 0;
}

//
// Print the a banner with time, date, and parallelization info (if any)
//
void
print_banner(char **argv,char *input_filename)
{
  time_t timenow = time(0);
  dft_log("\n******************************************************\n");
  dft_log("Current date and time: %s\n",ctime(&timenow));
#ifdef DFT_MPI
  dft_log("This is an MPI run:  N_Procs = %d\n",System::Get_N_Procs());
  dft_log("Setting NODE_ID dependent seed\n");
  srand(System::Get_procID()*1234+5678);//include time in it later
#endif
#ifdef DFT_THREAD
  dft_log("This is a thread run:  N_threads = %d\n",System::Get_N_threads());
#endif
  dft_log("%s: input file is '%s'\n",argv[0],input_filename);
}

//
// Print out help information (e.g. command line options)
//
void
help(char *name)
{
  // print out some help information
  dft_log("\nUsage: %s [options]\n",name);
  dft_log("\n\tPerforms ab-initio density-functional calculations.\n\n");
  dft_log("options:\n\n");
  dft_log("\t-h\t\thelp (this output)\n");
  dft_log("\t-i <filename>\tspecifies input file, default = 'dft.in'\n");
  dft_log("\t-o <filename>\tspecifies output log file, default = stdout\n");
#ifdef DFT_THREAD 
  dft_log("\t-n <integer>\tspecifies the number of threads, default = 1\n"); 
#endif 
  dft_log("\t-t\t\tprints an example input file to stdout and exits.\n");
  dft_log("\n");
}

//
// Parse command line.  Command line options currently
// are -i, -o, -t, and -h.
//
void
parse_command_line(int argc, char **argv, 
		   char *input_filename,
		   char *log_filename)
{
  // Logs to stdout.  We need this to report
  // basic errors and problems on the command line.
  // Make the system logging point to it for now.
  Output helplog;
  System::global_log = &helplog;

  int k=0;
  while(1)
    {
      // If we get to the end of the command line with no problems,
      // we are done!
      if (++k==argc)
	return;
      // produce a template input file.
      else if (MATCH(argv[k],"-t"))
	{
	  Everything tempeverything;
	  print_default_template(dft_global_log,tempeverything);
	  break;
	}
      // override default output file.
      else if (MATCH(argv[k],"-o"))
	{
	  if (++k<argc)
	    strcpy(log_filename,argv[k]);
	  else
	    {
	      dft_log("\nOption -o requires an argument.\n\n");
	      break;
	    }
	}
      // override default input file.
      else if (MATCH(argv[k],"-i"))
	{
	  if (++k<argc)
	    strcpy(input_filename,argv[k]);
	  else
	    {
	      dft_log("\nOption -i requires an argument.\n\n");
	      break;
	    }
	}
#ifdef DFT_THREAD
      // override the default number of threads.
      else if (MATCH(argv[k],"-n"))
	{
	  if (++k<argc)
	    System::N_threads=atoi(argv[k]);
	  else
	    {
	      dft_log("\nOption -n requires an integer argument.\n\n");
	      break;
	    }
	}
#endif
      // Help or known option
      else
	{
	  // If the option is known, complain
	  if (!MATCH(argv[k],"-h"))
	    dft_log("\nUnknown option: %s.\n",argv[k]);
	  break;
	}
    }

  // If we got here, print help and die
  help(argv[0]);
  System::GlobalFinalize();
  exit(0);
}


// ipd: what is the argument type? template this? noooooooooooooooooooo -csg
// typedef rulez!!!!! ivan is a silly boy
// ipd: yeah dude, it rulez, but it's k&r, wanna switch to fortran???

void
compare(GenericColumn &c1, GenericColumn &c2,
        char name1[20]="c1", char name2[20]="c2")
{
  dft_log("\nComparing:  %s (%s , %s)   ",
          name1,c1.getrepresentation(),c1.getembedding()); 
  dft_log("and    %s (%s , %s)\n",
          name2,c2.getrepresentation(),c2.getembedding()); 

  double m1, m2, m;

  m1=sqrt(abs(c1^c1));
  m2=sqrt(abs(c2^c2));
  {
    GenericColumn d(c1);
    d -= c2;
    m=sqrt(abs(d^d));
  }


  dft_log("|%s|\t= %le\n", name1, m1);
  dft_log("|%s|\t= %le\n", name2, m2);
  dft_log("|%s - %s|\t= %le\n", name1, name2, m);
  // this is nasty, maybe we only want to print it if m1 and m2 are not small?
  // but its not important...
  complex ipd = (c1^c2);
  complex nipd = ipd/(1.e-100+m1*m2);
  dft_log("|%s^%s|\t= %le + %le i\n[|%s^%s|]\t= %le + %le i\n\n",
          name1, name2, ipd.x, ipd.y, name1, name2, nipd.x, nipd.y);
}


void
optest(Everything &e)
{
  int i, j, q;
   
  dft_log("accessing via ColumnBundle.cols[i].data\n");
  for (q = 0; q < e.elecvars.nstates; q++)
    for(i = 0; i < e.elecvars.states[q].Y.my_ncols; i++)
      for(j = 0; j < e.elecvars.states[q].Y.col[i]->length; j++){
        (void)e.elecvars.states[q].Y.col[i]->data.d[j].x;
        (void)e.elecvars.states[q].Y.col[i]->data.d[j].y;
      }


  RealSpaceScalarFieldColumn &n = e.elecvars.n;
  RealSpaceScalarFieldColumn &d = e.elecvars.d;
  CoeffSpaceScalarFieldColumn m(d), mmm(d);
  RealSpaceScalarFieldColumn mm(n);

  real mag = 0.0;


  /**************************************
   * LJ1 = 0
   **************************************/
  dft_log("=============================\n");
  dft_log("LJ1 = 0 Test\n");
  dft_log("=============================\n");

  n.setmodones();

  dft_log("Calling J(n)\n"); dft_log_flush();
  apply_J(n,m);
  dft_log("Calling L(m)\n"); dft_log_flush();
  apply_L(m,mmm);

  compare(m,mmm, "n", "LJn");


  /**************************************
   * IJ = 1
   **************************************/
  dft_log("=============================\n");
  dft_log("IJ = 1 Test\n");
  dft_log("=============================\n");

  n.randomize();

  dft_log("Calling J(n)\n"); dft_log_flush();
  apply_J(n,m);
  dft_log("Calling I(m)\n"); dft_log_flush();
  apply_I(m,mm);

  compare(n, mm, "n", "IJn");

  /**************************************
   * JI = 1
   **************************************/
  dft_log("=============================\n");
  dft_log("JI = 1 Test\n");
  dft_log("=============================\n");


  d.randomize();

  dft_log("Calling I(d)\n"); dft_log_flush();
  apply_I(d,mm);
  dft_log("Calling J(mm)\n"); dft_log_flush();
  apply_J(mm,mmm);


  compare(d, mmm, "d", "JId");

  /**************************************
   * Idag = I^
   **************************************/
  dft_log("=============================\n");
  dft_log("Idag = I^  Test\n");
  dft_log("=============================\n");

  d.randomize();
  dft_log("d: %s %s\n", d.getrepresentation(), d.getembedding());
  m.randomize();
  dft_log("m: %s %s\n", m.getrepresentation(), m.getembedding());
  dft_log("Calling I(m)\n"); dft_log_flush();
  apply_I(m,m);// must be in real space
  dft_log("Calling I(d)\n"); dft_log_flush();
  apply_I(d,mm);
  dft_log("Calling Idag(m)\n"); dft_log_flush();
  apply_Idag(m,mmm);

  {
    complex diff=(m^mm)/(mmm^d)-1;
    dft_log("Idag is dag of I to %le + %le i\n\n",diff.x, diff.y );
  }

  
    
  /**************************************
   * Jdag = J^
   **************************************/
  dft_log("=============================\n");
  dft_log("Jdag = J^  Test\n");
  dft_log("=============================\n");

  d.randomize();
  dft_log("d: %s %s\n", d.getrepresentation(), d.getembedding());
  m.randomize();
  dft_log("m: %s %s\n", m.getrepresentation(), m.getembedding());
  dft_log("Calling J(m)\n"); dft_log_flush();
  apply_J(m,mmm);// must be in real space
  dft_log("Calling Jdag(d)\n"); dft_log_flush();
  apply_Jdag(d,mm);


  {
    complex diff=(m^mm)/(mmm^d)-1;
    dft_log("Jdag is dag of J to %le + %le i\n\n",diff.x, diff.y );
  }


  
  /**************************************
   * J1^OJ1 = Vol
   **************************************/
  dft_log("=============================\n");
  dft_log("J1^OJ1 = Vol Test\n");
  dft_log("=============================\n");

  for(i = 0; i < n.length; i++){
    n.data.d[i].x = 1;
    n.data.d[i].y = 0;
  }

  dft_log("Calling J(n)\n"); dft_log_flush();
  apply_J(n,m);
  dft_log("Calling O(m)\n"); dft_log_flush();
  apply_O(m,mmm);

  compare(m,mmm, "J1", "OJ1");
  



  /**************************************
   * J1^OJ1 = Vol on ColumnBundles
   **************************************/
  dft_log("==================================\n");
  dft_log("(J1)^OJ1 = Vol on ColumnBundle Test\n");
  dft_log("==================================\n");

  for(i = 0; i < n.length; i++){
    n.data.d[i].x = 1;
    n.data.d[i].y = 0;
  }

  mag = 0.0;
  for(i = 0; i < n.length; i++)
    mag += abs(n.data.d[i]);
  dft_log("\n mag[ n ] = %lf\n",mag);

  ColumnBundle &C = e.elecvars.states[0].C; 
  ColumnBundle &Y = e.elecvars.states[0].Y; 
  for(i = 0; i < C.my_ncols; i++){
    dft_log("Calling J(n)\n");
    apply_J(n, *(C.col[i]));
  }

  dft_log("Calling O(Y)\n"); dft_log_flush();
  O(C,Y);

  
  dft_log("Calling C^Y\n"); dft_log_flush();
  Matrix COC = C^Y;
  dft_log("C^OC : \n");
  COC.print();

  dft_log("\n");

  /**************************************
   * Orthonorm test (on ColumnBundles)
   **************************************/

  dft_log("==================================\n");
  dft_log("Orthonorm test (on ColumnBundles)\n");
  dft_log("==================================\n");

  C.randomize();
  COC=C^O(C);
  dft_log("\nOverlaps before:\n");
  
  COC.print();

  COC.hermetian = TRUE;

  Matrix Umh=Uminusonehalf(COC,e.elecvars.states[0].W,e.elecvars.states[0].mu);
  Y=C*Umh;


  dft_log("\nOverlaps after:\n");
  (Y^O(Y)).print();
  
  dft_log("\n");


  /**************************************
   * Test on B*M algebra: (B*B)**2=B*B^
   **************************************/

  dft_log("==================================\n");
  dft_log("B*M algebra test (B*B^)**2=B*B^  [idenpotency]\n");
  dft_log("==================================\n");

  //get idenpotent projector w/o O:
  COC = C^C;
  COC.hermetian = TRUE;
  Y = C*Uminusonehalf(COC,e.elecvars.states[0].W,e.elecvars.states[0].mu);

  dft_log("Projector idenpotency check:\n");
  (Y^Y).print();
  
  int col_num_diff = 3; // whatever to make it different
  ColumnBundle B(Y.col_length, Y.tot_ncols-col_num_diff, Y.col[0]->basis);
  // just put junk there, we're testing B*M
  B.randomize();

        
  ColumnBundle B1(B);
  B1 = Y*(Y^B);

  // project again
  B = Y*(Y^B1);

  dft_log("\n(PB)^(PB):\n");
  (B^B).print();

  //now compare
  B-=B1;

  dft_log("\nB^(PP-P)^(PP-P)B:\n");
  (B^B).printe();

  dft_log("\n");
}


void 
Etest(Everything &e)
{
//      calc_UVCn(e.elecinfo, e.elecvars, e.lattice, e.symm);

//      solve_poisson(e.elecvars.n,e.elecvars.d);
  calc_UVCn_d(e);
  
  //e.elecvars.n.read("n_new");
  //e.elecvars.n_up.read("n_up_new");
  //e.elecvars.n_dn.read("n_dn_new");

  calc_all_energies(e);

  //e.elecvars.n_up.write("n_up");
  //e.elecvars.n_dn.write("n_dn");

  e.energies.print(System::global_log);
}