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

  Class for canonical MC simulation in a restricted small cell (Fast)
  (Acceptance for multicomponent)

        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 <math.h>

#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"

#include "surf/cellBank.cpp"


typedef std::vector<std::vector<std::vector<std::vector<int> > > > vector4;
typedef std::map<std::vector<int>, double> mapVecD;

class GSMC{

protected:

  std::vector<Cluster> clusterArray;
  std::vector<Atom> atomArrayMC;
  std::vector<std::vector<int> > splitCellArray;
  vector4 coordination;

  std::vector<int> tempMC;
  int step1MC;
  int step2MC;

  /// MC Fast ///
  int numCell;
  int numBaseAtom;
  ///////////////

  // multi
  std::vector<std::vector<double> > basisFunc;
  
  std::ofstream cmcFile;

  // multi
  void setClusterArray(std::vector<Atom>& atomArrayCorre, std::vector<Atom>& atomArrayCorreLattice,
                       double& supercellSize, double & supercellSizeCorre,
                       std::map<int, mapVecD>& eci);
  void setAtomArray(std::vector<Atom>& atomArraySupercell, std::vector<Atom>& atomArrayCorreLattice,
		    const std::vector<int>& nCell);
  // multi
  void setPosition(int& natom_corre, const std::vector<std::vector<int> >& positionOut,
                   std::map<int, mapVecD>& eci);

  virtual void monteCarlo(const std::vector<int>& temp, int nStepAnn, int nStepEqv);
  // multi
  virtual mapVecD calcSigmaProd(std::vector<std::vector<int> > atomVector, int flipAtom, int flipSpin,
                                int anotherAtom, std::vector<std::vector<int> >& basisIn);

  // multi (new method)
  mapVecD multiplyMapVec(mapVecD mapVec, double scalar);
  void initializeMapVec(mapVecD& mapVec, std::vector<std::vector<int> >& vecs, double scalar);
  mapVecD addMapVec(mapVecD mapVec, mapVecD mapVecAdd);
  void selfAddMapVec(mapVecD& mapVec, mapVecD mapVecAdd);
  mapVecD diffMapVec(mapVecD mapVec, mapVecD mapVecAdd);
  mapVecD divideMapVec(mapVecD mapVec, mapVecD mapVec2, double scalar);
  
public: 
  //multi
  GSMC(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, mapVecD>& eci,
       const std::vector<std::vector<int> >& positionOut, std::vector<int>& splitCellIn);
  ~GSMC();

  virtual void start();


};

// multi
GSMC::GSMC(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, mapVecD>& eci,
           const std::vector<std::vector<int> >& positionOut, std::vector<int>& splitCellIn){

  clock_t startTime;
  clock_t currentTime;
  clock_t endTime;
  startTime = clock();

  MkSupercell mkcorre(atomArray, 2, 2, 2), mkcorre_lattice(atomArrayLattice, 2, 2, 2);
  std::vector<Atom> atomArrayCorre = mkcorre.getSupercell();
  std::vector<Atom> atomArrayCorreLattice = mkcorre_lattice.getSupercell();
  
  double sizeSupercell = atomArraySupercell.size() / atomArray.size();
  double sizeSupercellCorre = atomArraySupercell.size() / atomArrayCorre.size();
  
  std::cout << "setClusterArray ... " ;
  setClusterArray(atomArrayCorre, atomArrayCorreLattice, sizeSupercell, sizeSupercellCorre, eci);
  currentTime = clock();
  std::cout << "done.   Total elapsed time:" << (currentTime - startTime) / CLOCKS_PER_SEC << " sec" << std::endl;
  
  setAtomArray(atomArraySupercell, atomArrayCorreLattice, nCell);
  
  int natom_corre = atomArrayCorreLattice.size();
  
  
  /// MC Fast ///
  std::cout << "set Cell Bank ...";
  CellBank cellBank(atomArrayMC, splitCellIn);
  splitCellArray = cellBank.getCellBank();
  numCell = 1;
  std::vector<std::vector<int> >::iterator its;
  for (int i = 0; i < splitCellIn.size(); ++i){
    numCell *= splitCellIn[i];
  }
  numBaseAtom = atomArrayMC.size() / numCell;
  if (numBaseAtom != splitCellArray.size()){
    std::cout << "     Error: cannot divide cell. Please check the name: SPLITCELL." << std::endl;
    exit(0);
  }
  its = splitCellArray.begin();
  while (its != splitCellArray.end()){
    if ((*its).size() != numCell){
      std::cout << "     Error: cannot find equivalent atoms. Please check the name: SPLITCELL." << std::endl;
      exit(0);	    
    }
    ++its;
  }
  std::cout << "done." << std::endl;
  ///////////////
  
  
  std::cout << "setPosition ... " ;	
  setPosition(natom_corre, positionOut, eci);
  currentTime = clock();
  std::cout << "done.   Total elapsed time:" << (currentTime - startTime) / CLOCKS_PER_SEC << " sec" << std::endl;
  
  ///  ///
  step1MC = nStepAnn * splitCellArray.size();
  step2MC = nStepEqv * splitCellArray.size();
  for (int j = 0; j < temp.size(); ++j) tempMC.push_back(temp[j]);
  //////////
  
}

GSMC::~GSMC(){}

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

///// multi (new method) /////
std::map<std::vector<int>, double> GSMC::multiplyMapVec(mapVecD mapVec, double scalar){
  mapVecD results;
  mapVecD::iterator it;
  it = mapVec.begin();
  while (it != mapVec.end()){
    results[(*it).first] = (*it).second * scalar;
    ++it;
  }
  return results;
}
//////////////////////////////

// multi
void GSMC::setClusterArray(std::vector<Atom>& atomArrayCorre, std::vector<Atom>& atomArrayCorreLattice,
                           double& supercellSize, double & supercellSizeCorre,
                           std::map<int, mapVecD>& eci){
  
  ReadInput input;
  int numSublattice = input.getNumSublattice();
  const int maxCluster = input.getMaxCluster();
  const double tol = input.getTol();
  const double truncatePair = input.getTruncatePair();
  const double truncateTriangle = input.getTruncateTriangle();
  const double truncateTetrahedron = input.getTruncateTetrahedron();
  const double truncateMultiplet = input.getTruncateMultiplet();
  std::vector<std::string> nameArray = input.getAtomName();
  std::vector<int> spinArray = input.getSpinArray();
  int clusterOut = 1;
  
  int natom_index = atomArrayCorreLattice.size();
  EquivAtom equiv(atomArrayCorreLattice, natom_index);
  std::set<std::set<int> > equivAtom = equiv.getEquivAtom();
  std::vector<int> independentAtom = equiv.getIndependentAtom();
  std::vector<std::vector<int> > equivArray = equiv.getEquivArray();
  equivArray.resize(atomArrayCorreLattice.size());
  int unitCellAtom = atomArrayCorre.size()/8;
  
  
  MultiComponent multi(maxCluster, spinArray);
  std::vector<std::map<std::vector<int>, std::vector<double> > > correlationArray;
  GetCorrelation correlation(atomArrayCorre, atomArrayCorreLattice,
                             independentAtom, equivArray, truncatePair,truncateTriangle,
                             truncateTetrahedron, truncateMultiplet, tol, maxCluster,
                             unitCellAtom, clusterOut, equivAtom, multi);
  correlationArray = correlation.getCorrelation();
  // multi
  basisFunc = multi.getBasis();
  
  ///// multi /////
  std::map<int, mapVecD >::iterator ite = eci.begin();
  std::map<std::vector<int>, std::vector<double> >::iterator it2;
  std::map<std::vector<int>, double> sigmaProd, c0tmp, stmp, ttmp, utmp, wtmp;
  std::vector<int> vtmp, cvtmp;
  vtmp.push_back(0);
  stmp.insert(std::pair<std::vector<int>, double>(vtmp, 1.0));
  ttmp = multiplyMapVec((*ite).second, supercellSize);

  cvtmp.push_back(0);
  c0tmp[cvtmp] = 1;
  Cluster tmp( stmp , c0tmp,  0, ttmp );
  clusterArray.push_back(tmp);
  ++ite;
  
  while ( ite != eci.end()){
    int index = (*ite).first;
    if (index !=0){
      it2 = (correlationArray[index-1]).begin();
      std::map<std::vector<int>, double> mvtmp, nclusterTmp;
      double sigma, nCluster;
      while (it2 != (correlationArray[index-1]).end()){
        std::vector<int> basis = (*it2).first;
        sigma = ((*it2).second)[0];
        nCluster = ((*it2).second)[1];
        mvtmp[basis] = sigma;
        nclusterTmp[basis] = nCluster;
        ++it2;
      }
      ttmp = multiplyMapVec(mvtmp, supercellSizeCorre);
      utmp = multiplyMapVec((*ite).second, supercellSize);
      wtmp = multiplyMapVec(nclusterTmp, supercellSizeCorre);
      Cluster tmp( ttmp, wtmp, index, utmp );
      clusterArray.push_back(tmp);
      ++ite;
    }
  }
  /////////////////
}

void GSMC::setAtomArray(std::vector<Atom>& atomArraySupercell, std::vector<Atom>& atomArrayCorreLattice,
                        const std::vector<int>& nCell){
  
  std::vector<Atom> atomArrayTmp, atomArrayTmp2;
  atomArrayTmp = atomArrayCorreLattice;
  
  for (int i = 0; i < atomArraySupercell.size(); ++i){
    vector3D<double> positionSupercell = atomArraySupercell[i].getPosition();
    for (int j = 0; j < atomArrayCorreLattice.size(); ++j){
      vector3D<double> position = atomArrayCorreLattice[j].getPosition();
      double x = position.x * 2 / nCell[0];
      double y = position.y * 2 / nCell[1];
      double z = position.z * 2 / nCell[2];
      if (fabs(positionSupercell.x - x) < 1e-6
          and fabs(positionSupercell.y - y) < 1e-6
          and fabs(positionSupercell.z - z) < 1e-6){
        atomArrayTmp[j] = atomArraySupercell[i];
        break;
      }
      if (j == atomArrayCorreLattice.size() - 1)
        atomArrayTmp2.push_back(atomArraySupercell[i]);
    }
    
  }
  
  for (int i = 0; i < atomArrayTmp.size(); ++i) atomArrayMC.push_back(atomArrayTmp[i]);
  for (int i = 0; i < atomArrayTmp2.size(); ++i) atomArrayMC.push_back(atomArrayTmp2[i]);
  
}

// multi
void GSMC::setPosition(int& natom_corre, const std::vector<std::vector<int> >& positionOut,
                       std::map<int, mapVecD>& eci){ 
  
  EquivAtom equiv(atomArrayMC, natom_corre);
  std::set<std::set<int> > equivAtom = equiv.getEquivAtom();
  std::vector<int> independentAtom = equiv.getIndependentAtom();
  std::vector<std::vector<int> > equivArray = equiv.getEquivArray();
  equivArray.resize(natom_corre);
  
  std::vector<std::vector<int> > position;
  std::vector<int> pointIndex;
  
  int cutIndex = independentAtom.size();
  // multi
  std::map<int, mapVecD >::iterator it = eci.begin();
  while ( it != eci.end()){
    int index = (*it).first;
    if ( index > cutIndex ) position.push_back(positionOut[index - cutIndex - 1]);
    ++it;
    if ( index <= cutIndex and index != 0 ) pointIndex.push_back(index);
  }
  
  std::vector<std::set<std::set<int> > > clusterPosition;
  
  std::set<std::set<int> >::iterator it0;
  std::set<int>::iterator it00;
  std::set<std::set<int> > tmp0;
  std::set<int> tmp00;
  ///  ///
  for (int i = 0; i < pointIndex.size(); ++i){
    //////////
    it0 = equivAtom.begin();
    for (int j = 1; j < pointIndex[i]; ++j) ++it0;
    it00 = (*it0).begin();
    while( it00 != (*it0).end()){
      tmp00.insert(*it00);
      tmp0.insert(tmp00);
      tmp00.clear();
      ++it00;
    }
    clusterPosition.push_back(tmp0);
    tmp0.clear();
  }
  
  for (int i = 0; i < position.size(); ++i){
    for (int j = 0; j < equivArray[0].size(); ++j){
      for (int k = 0; k < position[i].size(); ++k){
        tmp00.insert(equivArray[position[i][k]][j]);
      }
      tmp0.insert(tmp00);
      tmp00.clear();
    }
    clusterPosition.push_back(tmp0);
    tmp0.clear();
  }
  equivArray.clear();
  
  coordination.resize(atomArrayMC.size());
  for (int i = 0; i < coordination.size(); ++i){
    coordination[i].resize(clusterPosition.size());
  }
  
  for (int i = 0; i < clusterPosition.size(); ++i){
    it0 = clusterPosition[i].begin();
    while ( it0 != clusterPosition[i].end()){
      std::vector<int> tmpVec;
      it00 = (*it0).begin();
      while ( it00 != (*it0).end()){
        tmpVec.push_back(*it00);
        ++it00;
      }
      
      it00 = (*it0).begin();
      while ( it00 != (*it0).end()){
        coordination[*it00][i].push_back(tmpVec);
        ++it00;
      }
      ++it0;
      tmpVec.clear();
    }
  }
  
}


void GSMC::monteCarlo(const std::vector<int>& temp, int nStepAnn, int nStepEqv){
  
  const double k_B = 0.86171 * pow (10,-4);
  unsigned long init = time(NULL) + getpid();
  init_genrand(init);
  
  // multi
  std::vector<mapVecD> correlationArray;
  correlationArray.resize(clusterArray.size());
  
  
  long double energy(0), energyAve(0), energySquare(0);
  int natom = splitCellArray.size();

  ///// multi /////
  std::cout << std::endl << "Cluster, [basisIndex, correlation, ECI]: " << std::endl;
  for (int i = 0; i < clusterArray.size(); ++i) {
    int cIndex = clusterArray[i].getIndex();
    mapVecD nClusterMap = clusterArray[i].getNumCluster();
    mapVecD sigmatmp = clusterArray[i].getSpin();
    mapVecD::iterator itsigma = sigmatmp.begin();
    std::cout << cIndex << "    ";
    while (itsigma != sigmatmp.end()){
      std::vector<int> basis;
      basis = (*itsigma).first;
      double sigmaProd = (*itsigma).second;
      double eciTmp = (clusterArray[i].getECI())[basis];
      double nCluster = nClusterMap[basis];
      energy += eciTmp * sigmaProd / nCluster;
      for (int i = 0; i < basis.size(); ++i){
        std::cout << basis[i];
      }
      std::cout << "  " << sigmaProd / nCluster << ", " << eciTmp << "  ";
      ++itsigma;
    }
    std::cout << std::endl;
  }
  /////////////////
  
  std::cout << "E (initial) = "<< energy << std::endl;

  // BEGIN: COMMENT OUT 

    
  // 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){
      
      int atom1, atom2;
      int spin1, spin2;
     
      /// Select 2 atoms in the "base" cell randomly
      atom1 = genrand_int32() % natom;
      atom2 = genrand_int32() % natom;
      spin1 = atomArrayMC[splitCellArray[atom1][0]].getSpin();
      spin2 = atomArrayMC[splitCellArray[atom2][0]].getSpin();
      
      if ( spin1 != spin2 )  { 

        /// multi
        mapVecD spinOld, spinNew, diffSpin;
        std::vector<mapVecD> diffSpinArray;
        long double diffEnergy = 0;
        
        /// Loop for respective clusters
        for (int i = 0; i < clusterArray.size() - 1; ++i){
          
          /// Flipping spins for ALL the cells successively ///
          mapVecD diffSpinSum;
          int atomC1, atomC2, spinC1, spinC2;
          // multi
          std::vector<std::vector<int> > basisIndex;
          basisIndex = clusterArray[i+1].getBasisIndex(); // Because i=0 is empty cluster, i=0 should be omitted
                                                          // when taking energy difference.
          initializeMapVec(diffSpinSum, basisIndex, 0.0);
          
          for (int j = 0; j < numCell; ++j){
            atomC1 = splitCellArray[atom1][j];
            atomC2 = splitCellArray[atom2][j];
            spinC1 = atomArrayMC[atomC1].getSpin();
            spinC2 = atomArrayMC[atomC2].getSpin();

            mapVecD spTmp1, spTmp2, spTmp3;

            // Coordination begin from point cluster (not empty). So, array index is [i].
            spTmp1  = calcSigmaProd(coordination[atomC1][i], atomC1, spinC1, atomC2, basisIndex);
            spTmp2  = calcSigmaProd(coordination[atomC2][i], atomC2, spinC2, atomC1, basisIndex);
            spinOld = addMapVec(spTmp1, spTmp2);

            spTmp1  = calcSigmaProd(coordination[atomC1][i], atomC1, spinC2, atomC2, basisIndex);
            spTmp2  = calcSigmaProd(coordination[atomC2][i], atomC2, spinC1, atomC1, basisIndex);
            spinNew = addMapVec(spTmp1, spTmp2);
            
            spTmp3 = diffMapVec(spinNew, spinOld);
            
            diffSpinSum = addMapVec(diffSpinSum, spTmp3);
            
            /// At this time, spins of 2 atoms should be flipped. ///
            atomArrayMC[atomC1].setSpin(spinC2);
            atomArrayMC[atomC2].setSpin(spinC1);
            //////////
            
          }
          
          /// Putting back the spins to original value ///
          for (int k = 0; k < numCell; ++k){
            atomC1 = splitCellArray[atom1][k];
            atomC2 = splitCellArray[atom2][k];
            spinC1 = atomArrayMC[atomC1].getSpin();
            spinC2 = atomArrayMC[atomC2].getSpin();
            atomArrayMC[atomC1].setSpin(spinC2);
            atomArrayMC[atomC2].setSpin(spinC1);
          }
          ///
          
          diffSpinArray.push_back(diffSpinSum);

          
          /// calc diff energy ///
          mapVecD nClusterMap2 = clusterArray[i+1].getNumCluster();
          mapVecD::iterator itsigma2 = diffSpinSum.begin();
          double energy2 = 0.0;
          while (itsigma2 != diffSpinSum.end()){
            std::vector<int> basis2;
            basis2 = (*itsigma2).first;
            double sigmaProd2 = (*itsigma2).second;
            double eciTmp2 = (clusterArray[i+1].getECI())[basis2];
            double nCluster2 = nClusterMap2[basis2];
            energy2 += eciTmp2 * sigmaProd2 / nCluster2;
            ++itsigma2;
          }
          
          diffEnergy +=  energy2;
          ///////////////////////
          
          
        }
        
        double tmpProb = -diffEnergy / k_B / tempTrial;
        double probability = exp ( tmpProb );
        double cut = genrand_real2();
        
        int atomC1, atomC2, spinC1, spinC2;
        if (probability > cut ) {
          energy += diffEnergy;
          /// Flipping All the atoms in each cells
          for (int k = 0; k < numCell; ++k){
            atomC1 = splitCellArray[atom1][k];
            atomC2 = splitCellArray[atom2][k];
            spinC1 = atomArrayMC[atomC1].getSpin();
            spinC2 = atomArrayMC[atomC2].getSpin();
            atomArrayMC[atomC1].setSpin(spinC2);
            atomArrayMC[atomC2].setSpin(spinC1);
          }
          ///

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

        if (iterTemp == temp.size()){
          for (int i = 1; i < clusterArray.size(); ++i){
            mapVecD spinN = clusterArray[i].getSpin();
            mapVecD numClusterN = clusterArray[i].getNumCluster();
            mapVecD mv1, mv2;
            mv1 = divideMapVec(spinN, numClusterN, nStepEqv);
            selfAddMapVec(correlationArray[i], mv1);
            //////////
            
          }

        }
        
        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("gsmc.out", std::ios::out);
  cmcFile << "------------------------------------------" << std::endl;
  cmcFile << "	Averages of correlation functions     " << std::endl;
  cmcFile << "------------------------------------------" << std::endl;
  energyAve = 0;
  
  ///// multi /////////////////////////////////////////////
  mapVecD cctmp;
  std::vector<int> vvtmp;
  vvtmp.push_back(0);
  cctmp[vvtmp] = 1;
  correlationArray[0] = cctmp;

  for (int i = 0; i < clusterArray.size(); ++i) {
    mapVecD sigmatmp = correlationArray[i];
    mapVecD::iterator itsigma = sigmatmp.begin();
    int index = clusterArray[i].getIndex();
    cmcFile << index << " : ";
    while (itsigma != sigmatmp.end()){
      std::vector<int> basis;
      basis = (*itsigma).first;
      double sigmaProd = (*itsigma).second;
      double eciTmp = (clusterArray[i].getECI())[basis];
      energyAve += eciTmp * sigmaProd;
      for (int i = 0; i < basis.size(); ++i){
        cmcFile << basis[i];
      }
      cmcFile << "  " << sigmaProd  << "    ";
      ++itsigma;
    }
    cmcFile << std::endl;
  }

  ///////////////////////////////////////////////





  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();


   //
   // END: COMMENTOUT
    
}


mapVecD GSMC::calcSigmaProd(std::vector<std::vector<int> > atomVector,
                           int flipAtom, int flipSpin, int anotherAtom,
                           std::vector<std::vector<int> >& basisIndex){

  mapVecD sigmaProdTotal;
  initializeMapVec(sigmaProdTotal, basisIndex, 0.0);
  
  // Loop for each clusters (same type)
  for (int i = 0; i < atomVector.size(); ++i){
    mapVecD sigmaProd;

    double spin;
    // Loop for each basis index set: exp -> basisIndexSet = [[1, 1, 1], [1, 1, 2], ..] 
    for (int p=0; p < basisIndex.size(); ++p){
      std::vector<int> basis = basisIndex[p];
      double sigma;
      double sigmaSum = 0.0;
      sigma = 1.0;
      // Size of basis index and order of cluster must match.
      // (Otherwise, this loop causes segmentation fault due to null-pointer exception)
      // Loop for basis: exp -> basis = [1, 1, 2]
      for (int q=0; q < basis.size(); ++q){
        if (flipAtom == atomVector[i][q]){
          spin = flipSpin;
        }else if(anotherAtom == atomVector[i][q]){
          // initializeMapVec(sigmaProd, basisIndex, 0.0);
          sigma = 0.0;
          break;
        }else{
          spin = atomArrayMC[atomVector[i][q]].getSpin();
        }
        std::vector<double> Psi = basisFunc[basis[q]];
        double PsiValue = 0.0;
        for (int r=0; r < Psi.size(); ++r){
          PsiValue += Psi[r] * pow(spin, r);
        }
        sigma *= PsiValue;
      }
      sigmaSum += sigma;

      // next index after sorting: exp -> basis(orig) = [1, 1, 2] --> basis(next) = [1, 2, 1]
      while (next_permutation(basis.begin(), basis.end())){
        sigma = 1.0;
        for (int q=0; q < basis.size(); ++q){
          if (flipAtom == atomVector[i][q]){
            spin = flipSpin;
          }else if(anotherAtom == atomVector[i][q]){
            // initializeMapVec(sigmaProd, basisIndex, 0.0);
            sigma = 0.0;
            break;
          }else{
            spin = atomArrayMC[atomVector[i][q]].getSpin();
          }
          std::vector<double> Psi = basisFunc[basis[q]];
          double PsiValue = 0.0;
          for (int r=0; r < Psi.size(); ++r){
            PsiValue += Psi[r] * pow(spin, r);
          }
          sigma *= PsiValue;
        }
        sigmaSum += sigma;        
      }
      
      sort(basis.begin(), basis.end());
      sigmaProd[basis] = sigmaSum;
    }
    // Add sigmaProd(for all basis index, but for one cluster) to sigmaProdTotal(for all cluster(same type))
    sigmaProdTotal = addMapVec(sigmaProdTotal, sigmaProd);
  }
  return sigmaProdTotal;
}




mapVecD GSMC::addMapVec(mapVecD mapVec, mapVecD mapVecAdd){
  mapVecD mtmp;
  mapVecD::iterator it;
  it = mapVecAdd.begin();
  while (it != mapVecAdd.end()){
    std::vector<int> key = (*it).first;
    double value = (*it).second;
    double value2 = mapVec[key];
    mtmp[key] = value + value2;
    ++it;
  }
  return mtmp;
}

void GSMC::selfAddMapVec(mapVecD& mapVec, mapVecD mapVecAdd){
  mapVecD::iterator it;
  it = mapVecAdd.begin();
  while (it != mapVecAdd.end()){
    std::vector<int> key = (*it).first;
    double value2 = mapVecAdd[key];
    mapVecD::iterator it2 = mapVec.find(key);
    if (it2 == mapVec.end()){
      mapVec[key] = 0.0;
    }
    mapVec[key] += value2;
    ++it;
  }
}

mapVecD GSMC::diffMapVec(mapVecD mapVec, mapVecD mapVec2){
  mapVecD mtmp;
  mapVecD::iterator it;
  it = mapVec2.begin();
  while (it != mapVec2.end()){
    std::vector<int> key = (*it).first;
    double value2 = (*it).second;
    double value = mapVec[key];
    mtmp[key] = value - value2;
    ++it;
  }
  return mtmp;
}

mapVecD GSMC::divideMapVec(mapVecD mapVec, mapVecD mapVec2, double scalar){
  mapVecD mtmp;
  mapVecD::iterator it;
  it = mapVec2.begin();
  while (it != mapVec2.end()){
    std::vector<int> key = (*it).first;
    double value2 = (*it).second;
    double value = mapVec[key];
    mtmp[key] = value / value2 / scalar;
    ++it;
  }
  return mtmp;
}




void GSMC::initializeMapVec(mapVecD& mapVec, std::vector<std::vector<int> >& vec, double scalar){
  std::vector<std::vector<int> >::iterator it;
  it = vec.begin();
  while (it != vec.end()){
    mapVec[(*it)] = scalar;
    ++it;
  }
}