/******************************************************************************

  Class for grand canonical surface MC simulation 

 	Written by k-fleak

 ********************************************************************************/

#pragma warning( disable : 4786 )
#include <iostream>
#include <fstream>
#include <vector>
#include <set>
#include <map>
#include <time.h>
#include <unistd.h>

#include "cmcSurf.cpp"

/*
#include "mt/mt19937ar.h"
#include "mt/mt19937ar.c"

#include "vector3D.hpp"
#include "parsePoscar.cpp"
#include "atom.cpp"
#include "mkSupercell.cpp"
#include "readInput.cpp"
#include "cluster.cpp"
#include "equivAtom.cpp"
#include "multiComponent.cpp"
#include "getCorrelation.cpp"
*/

typedef std::vector<std::vector<std::vector<std::vector<int> > > > vector4;

class GCMCSURF : public CMCSURF{
private:
  //   Note: Match the sequence of element name in POTCAR CE.in, and MC.in (mu).
  std::map<int, double> muArray; // [spin, mu]
  std::map<int, int> numAtomArray; // [spin, numAtom]
  std::vector<int> spinArray;
  void monteCarlo(const std::vector<int>& temp, int nStepAnn, int nStepEqv);
  double calcSpin(std::vector<std::vector<int> > atomVector, int flipAtom, int flipSpin);

public: 
  GCMCSURF(std::vector<Atom>& atomArray, std::vector<Atom>& atomArrayLattice, 
       std::vector<Atom>& atomArraySupercell, const std::vector<int>& temp, const int& nStepAnn, 
       const int& nStepEqv, const std::vector<int>& nCell, std::map<int, double>& eci,
       const std::vector<std::vector<int> >& positionOut, 
       std::vector<int> numAtomsIn, std::vector<double> muIn, std::vector<double> MovAtomZPosIn) ;
  ~GCMCSURF();
  void start();
};


GCMCSURF::GCMCSURF(std::vector<Atom>& atomArray, std::vector<Atom>& atomArrayLattice, 
	   std::vector<Atom>& atomArraySupercell, const std::vector<int>& temp, const int& nStepAnn, 
	   const int& nStepEqv, const std::vector<int>& nCell, std::map<int, double>& eci,
	   const std::vector<std::vector<int> >& positionOut, 
	   std::vector<int> numAtomsIn, std::vector<double> muIn, std::vector<double> MovAtomZPosIn) 
  : CMCSURF(atomArray, atomArrayLattice, atomArraySupercell, temp, nStepAnn, 
	nStepEqv,  nCell, eci, positionOut, MovAtomZPosIn) {
  ReadInput input;
  spinArray = input.getSpinArray();

  std::cout << "gcmc constcutror" << std::endl;

  
  for ( int j = 0; j < muIn.size(); ++j) {
    muArray.insert(std::pair<int, double>(spinArray[j], muIn[j]));
    std::cout << "spin, mu: " << spinArray[j] << " " << muIn[j] << std::endl;
  }

  int cellSize, numAtomCell = 0;
  for ( int n = 0; n < numAtomsIn.size(); ++n) numAtomCell += numAtomsIn[n];
  cellSize = atomArrayMC.size() / numAtomCell;
  for ( int i = 0; i < numAtomsIn.size(); ++i) numAtomArray.insert(std::pair<int, int>(spinArray[i], numAtomsIn[i] * cellSize));
}


GCMCSURF::~GCMCSURF(){
}

void GCMCSURF::start(){
  clock_t startTime;
  clock_t endTime;
  startTime = clock();
  std::cout << "GC monteCarlo ... " ;
  monteCarlo(tempMC, step1MC, step2MC);
  endTime = clock();
  std::cout << "done.   GC monteCarlo elapsed time: " << (endTime-startTime) / CLOCKS_PER_SEC << " sec" << std::endl;  
}

void GCMCSURF::monteCarlo(const std::vector<int>& temp, int nStepAnn, int nStepEqv){
  
  std::cout << "GCMC" << std::endl;

  const double k_B = 0.86171 * pow (10,-4);
  unsigned long init = time(NULL) + getpid();
  init_genrand(init);

  std::map<int, double> numAtomAve; // [spin, numAtom]
  for ( int i = 0; i < spinArray.size(); ++i ) numAtomAve[spinArray[i]] = 0.0;

  std::vector<double> correlationArray;
  correlationArray.resize(clusterArray.size());
	
  long double energy(0), energyAve(0), energySquare(0);
  int natom = atomArrayMC.size();

  std::cout << "Initial correlation function: " << std::endl;
  for (int i = 0; i < clusterArray.size(); ++i) {
    energy += clusterArray[i].getECI() * clusterArray[i].getSpin() / clusterArray[i].getNumCluster();
    std::cout << clusterArray[i].getSpin() / clusterArray[i].getNumCluster() << std::endl;
  }

  std::cout << "E (initial) = "<< energy << std::endl;

  // BEGIN MAIN LOOP

  for (int iterTemp = 0; iterTemp < temp.size() + 1; ++iterTemp){
    int tempTrial, nStep;
    if (iterTemp == temp.size()){
      tempTrial = temp[temp.size()-1];
      nStep = nStepEqv;
    }
    else {
      tempTrial = temp[iterTemp];
      nStep = nStepAnn;
    }
	

    std::cout << "TEMPERATURE = "<< tempTrial << " " << "NSTEP = "<< nStep << std::endl;
    int stepSum = 1;

    while (stepSum <= nStep){

      // atom1: system  atom2: bath
      int atom1, atom2;
      int spin1, spin2;
      vector3D<double> pos1;
      int flag1;
      atom1 = genrand_int32() % natom;
      atom2 = genrand_int32() % spinArray.size();
      spin1 = atomArrayMC[atom1].getSpin();
      spin2 = spinArray[atom2];


      //      std::cout << "atom1,2, spin1,2 :" << atom1 << " " << atom2 << " " << spin1 << " " << spin2 <<  std::endl;

      /// surf ///
      pos1 = atomArrayMC[atom1].getPosition();
      flag1 = 0;
      if ( spin1 != spin2 ){
	for (int i = 0; i < slab1ZPos.size(); ++i ) {
	  if ( (pos1.z >= slab1ZPos[i] - eps) and (pos1.z <= slab1ZPos[i] + eps)) {
	    flag1 = 1;
	  } 
	}
	for (int i = 0; i < slab2ZPos.size(); ++i ) {
	  if ( (pos1.z >= slab2ZPos[i] - eps) and (pos1.z <= slab2ZPos[i] + eps)) {
	    flag1 = 1;
	  }
	}
      }
      ////////////

      if ( flag1 == 1) {
	double spinOld, spinNew, diffSpin;
	std::vector<double> diffSpinArray;
	long double diffEnergy = 0;

	for (int i = 0; i < clusterArray.size() - 1; ++i){
				
	  spinOld = calcSpin(coordination[atom1][i], atom1, spin1);
	  spinNew = calcSpin(coordination[atom1][i], atom1, spin2);
	  diffSpin = spinNew - spinOld;
				
	  diffSpinArray.push_back(diffSpin);
	  diffEnergy +=  clusterArray[i+1].getECI() * diffSpin / clusterArray[i+1].getNumCluster();

	}

	double tmpProb = - ( diffEnergy + muArray[spin1] - muArray[spin2] )  / k_B / tempTrial;
	double probability = static_cast<double>( numAtomArray[spin1] / (numAtomArray[spin2] + 1.0) )  * exp ( tmpProb );
	double cut = genrand_real2();


	//	std::cout << "diffE, dmu, tmpProb, prob, cut :" << diffEnergy << " " << muArray[spin1] - muArray[spin2] <<" "<< tmpProb << " " << probability << " " << cut << std::endl;
	//	std::cout << "numAtomArray[spin1], numAtomArray[spin2]: " << numAtomArray[spin1] << " " << numAtomArray[spin2] << std::endl;

	if (probability > cut ) {
	  energy += diffEnergy;
	  atomArrayMC[atom1].setSpin(spin2);
	  numAtomArray[spin1] -= 1;
	  numAtomArray[spin2] += 1;
	  

	  //	  std::cout << "Adopt. numAtomArray[spin1,2]: " << numAtomArray[spin1] << " " << numAtomArray[spin2] << std::endl;


	  for (int i = 0; i < clusterArray.size() - 1; ++i) clusterArray[i+1].addSpin(diffSpinArray[i]);
	}


	
	/*********************************************************************************/
	// check the correlation functions 
	/*********************************************************************************/
	
	//	for (int i = 0; i < clusterArray.size(); ++i){
	//		double spin_tmp = clusterArray[i].getSpin();
	//		double clusterNum_tmp = clusterArray[i].getNumCluster();
	//		if (fabs(spin_tmp) > clusterNum_tmp ){	
	//			std::cerr << "Error: Correlation Functions are over one!\n";
	//			std::cout << "Error: " << clusterArray[i].getSerialIndex() << " Correlation Functions are over one! " << std::endl;
	//	                exit(8);
	//		}
	//	}	
	/*********************************************************************************/
	// Average 
	/*********************************************************************************/

	if (iterTemp == temp.size()){
	  for (int i = 1; i < clusterArray.size(); ++i){
	    double spin = clusterArray[i].getSpin();
	    double numCluster = clusterArray[i].getNumCluster();
	    correlationArray[i] += spin / numCluster / nStepEqv;
	  }
	  for ( int i = 0; i < spinArray.size(); ++i ){
	    numAtomAve[spinArray[i]] += numAtomArray[spinArray[i]] / static_cast<double>(nStepEqv) ;
	  }
	}
	  diffSpinArray.clear();
	  ++stepSum;
	

      }
	
    } 
      std::cout << "E = "<< energy << " " << std::endl;
  }

  // END OF MAIN LOOP	
  std::cout << " End of loop " << std::endl;

  // WRITING OUTPUT FILE FOR cmc.out	

  cmcFile.open("gcmc.out",std::ios::out);
  //
  //	energyAve /= NstepEq;
  //	energySquare /= NstepEq;
  //	
  cmcFile << "------------------------------------------" << std::endl;
  cmcFile << "	Averages of correlation functions     " << std::endl;
  cmcFile << "------------------------------------------" << std::endl;
  energyAve = 0;
  correlationArray[0] = 1;
  for (int i = 0; i < clusterArray.size(); ++i){
    int index = clusterArray[i].getIndex();
    double numCluster = clusterArray[i].getNumCluster();
    //double correlation = correlationArray[i] / numCluster / NstepEq;
    double correlation = correlationArray[i];
    energyAve += clusterArray[i].getECI() * correlation;

    cmcFile << index << " : " << correlation << std::endl;
  }


  cmcFile << "------------------------------------------" << std::endl;
  cmcFile << "	        Averages of composition       " << std::endl;
  cmcFile << "------------------------------------------" << std::endl;
  std::map<int, double>::iterator it;
  it = numAtomAve.begin();
  while ( it != numAtomAve.end()){
    cmcFile << "spin: " << (*it).first << "   composition: " << (*it).second / atomArrayMC.size() << std::endl;
    ++it;
  }



  //
  cmcFile << "------------------------------------------" << std::endl;
  cmcFile << "	       Averages of energy             " << std::endl;
  cmcFile << "------------------------------------------" << std::endl;
  //	
  cmcFile.precision(10);
  //	cmcFile << "Number of steps of average = " << stepAverage << std::endl;
  cmcFile << "Average of energy = "<< energyAve << " eV" << std::endl;
  //	cmcFile << "Average of square of energy = "<< energySquare << " eV^2" << std::endl;
  //	//cmcFile << "Boltzmann constant = "<< k_B << std::endl;
  //	cmcFile << "Temperature = "<< tempEnd <nn< std::endl;
  //	long double squareTemp = k_B * tempEnd * tempEnd;
  //	long double squareEnergyAve = energyAve * energyAve;
  //	long double diffEnergyAve = energySquare - squareEnergyAve;
  //	cmcFile << "Specific heat at constant volume  = "<< diffEnergyAve / squareTemp  << std::endl;
  //
  cmcFile << std::endl;
  //
  cmcFile << "------------------------------------------" << std::endl;
  cmcFile << "	            Position                  " << std::endl;
  cmcFile << "------------------------------------------" << std::endl;
  //
  for (int i = 0; i < atomArrayMC.size(); ++i){
    vector3D<double> position = atomArrayMC[i].getPosition();
    cmcFile << position.x << " " << position.y << " " << position.z << " ";
    cmcFile << atomArrayMC[i].getSpin() << std::endl;

  }
  cmcFile << std::endl;
  cmcFile.close();
}

double GCMCSURF::calcSpin(std::vector<std::vector<int> > atomVector, int flipAtom, int flipSpin){

  double spinSum = 0;
  for (int i = 0; i < atomVector.size(); ++i){
    double spin = 1;
    for (int j = 0; j < atomVector[i].size(); ++j){
      if (flipAtom == atomVector[i][j]){
	spin *= flipSpin;
      }else{
	spin *= atomArrayMC[atomVector[i][j]].getSpin();
      }
    }
    spinSum += spin; 
  }
  return spinSum;
}