Manybody_Perturbative_Correction_Template.h

#ifndef QNANO_NEW_MANYBODY_PERTURBATIVE_CORRECTION_TEMPLATE_DEFINED_H
#define QNANO_NEW_MANYBODY_PERTURBATIVE_CORRECTION_TEMPLATE_DEFINED_H


template <class Tsolver> void Communicate_Complex_Vector_Template(Tsolver *solver, std::vector<std::complex<double> > &EV2){
}

#ifdef QNANO_NEW_SLEPC_SOLVER_DEFINED_H_
template <> void Communicate_Complex_Vector_Template<SLEPc_Solver>(SLEPc_Solver* solver, std::vector<std::complex<double> > &EV2){
  MPI_Allreduce(MPI_IN_PLACE, &EV2[0], EV2.size()*2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
}
#endif



template <class Tsolver> void Communicate_Append_Int_Vector_Template(Tsolver *solver, std::vector<int> &newstates_indices){
}
#ifdef QNANO_NEW_SLEPC_SOLVER_DEFINED_H_
template <> void Communicate_Append_Int_Vector_Template<SLEPc_Solver>(SLEPc_Solver *solver, std::vector<int> &newstates_indices){
  int MPI_rank, MPI_size;
  MPI_Comm_size(MPI_COMM_WORLD, &MPI_size);
  MPI_Comm_rank(MPI_COMM_WORLD, &MPI_rank);

  if(MPI_rank==0){
    for(size_t j=1; (int)j<MPI_size; j++){
      int vectorsize;
      MPI_Recv(&vectorsize, 1, MPI_INT, j, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      size_t oldsize=newstates_indices.size();
      newstates_indices.resize(oldsize+vectorsize);
      MPI_Recv(&newstates_indices[oldsize], vectorsize, MPI_INT, j, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    }
  }else{
    int vectorsize=newstates_indices.size();
    MPI_Send(&vectorsize, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
    MPI_Send(&newstates_indices[0], vectorsize, MPI_INT, 0, 0, MPI_COMM_WORLD);
  }
}
#endif

template <class Tsolver> void Calculate_Manybody_Perturbation(Parameter_Map &param_map, Tsolver *solver, Eigenvector_Printer *EVprinter, Manybody_Hamiltonian *hamil){

  std::string PT_outfile=param_map.get_as_string("PT_outfile","");
  if(PT_outfile=="")return;

  int nr_print=EVprinter->get_actual_nr_print();
  nr_print=param_map.get_as_double("nr_PT",nr_print);
  if(nr_print>EVprinter->get_n_converged())nr_print=EVprinter->get_n_converged();
  if(nr_print<1)return;

  bool should_print=solver->should_print();  //true for Eigen, or if MPI_rank=0 for SLEPc

  double small_denominator_threshold=param_map.get_as_double("small_denominator_threshold", 1e-5);


  //Hamiltonian in extended Hilbert space:
  if(should_print)std::cout<<"Setting up manybody Hamiltonian for perturbative corrections..."<<std::endl;

  int nr_e2=param_map.get_as_sizet("nr_e_PT",0);
  int nr_h2=param_map.get_as_sizet("nr_h_PT",0);
  if(nr_e2<hamil->hilbert->get_nr_e()){
    std::cerr<<"nr_e_PT < nr_e !"<<std::endl;
    exit(1);
  }
  if(nr_h2<hamil->hilbert->get_nr_h()){
    std::cerr<<"nr_h_PT < nr_h !"<<std::endl;
    exit(1);
  }
  int s_nr_e=param_map.get_as_sizet("selected_nr_e");
  if(!param_map.is_specified("selected_nr_e")){
    s_nr_e=hamil->hilbert->get_nr_e();
  }
  int s_nr_h=param_map.get_as_sizet("selected_nr_h");
  if(!param_map.is_specified("selected_nr_h")){
    s_nr_h=hamil->hilbert->get_nr_h();
  }


  Parameter_Map param_map2(param_map);
  param_map2.override_parameter("nr_e", Reader::toString(nr_e2));
  param_map2.override_parameter("nr_h", Reader::toString(nr_h2));
  param_map2.override_parameter("selected_nr_e", Reader::toString(s_nr_e));
  param_map2.override_parameter("selected_nr_h", Reader::toString(s_nr_h));

  Manybody_Hamiltonian hamil2(param_map2);
  {
    Selected_Manybody_Hilbert_Space hilbert2(*hamil->hilbert);
    hilbert2.nr_e=nr_e2;
    hilbert2.nr_h=nr_h2;
    hilbert2.add_singles();
    hilbert2.add_singles();

    hamil2.replace_hilbert_space(hilbert2);
  }

  
  if(should_print){
    std::cout<<"Dimension of original Hilbert space: "<<hamil->get_DIM()<<std::endl;
    std::cout<<"Dimension of new Hilbert space: "<<hamil2.get_DIM()<<std::endl;
  }
 
 
  //print new states if contribution is above a certail threshold
  Selected_Manybody_States_with_Hash_Table newstates; 
  std::vector<int> newstates_indices;
  double threshold=param_map.get_as_double("threshold");

  if(param_map.is_specified("new_states_file") &&
     param_map.is_specified("add_new_states_file") ){
    std::cerr<<"Please set either 'new_states_file' or 'add_new_states_file', but not both!"<<std::endl;
    exit(1);
  }

  int existing_states_in_new_states_file=0;
  std::string new_states_file=param_map.get_as_string("new_states_file");
  if(param_map.is_specified("add_new_states_file")){
    new_states_file=param_map.get_as_string("add_new_states_file");
    newstates.add_from_file(new_states_file, \
      hamil2.hilbert->get_rank_e(), hamil2.hilbert->get_rank_h(), 
      hamil2.hilbert->get_nr_e(), hamil2.hilbert->get_nr_h());
    existing_states_in_new_states_file=newstates.size();
  }

  bool print_new_states=false;
  if(new_states_file!=""){
    print_new_states=true;
  }


  // Prepare output stream
  std::ofstream ofs;
  if(should_print){
    ofs.open(PT_outfile.c_str());
    ofs<<"#DIM "<<hamil->get_DIM()<<" -> "<<hamil2.get_DIM()<<std::endl;
    std::cout<<std::setprecision(10);
    ofs<<std::setprecision(10);
  }
  

  Distribute_Block_Local *db=solver->new_Distribute_Block_Local(hamil->get_DIM());
  Distribute_Block_Local *db2=solver->new_Distribute_Block_Local(hamil2.get_DIM());
  Communicator *comm=solver->new_Communicator();
 

  //Handle EVs block-wise to save memory
  size_t nr_blocks=param_map.get_as_sizet("nr_blocks",1);
  if(param_map.is_specified("EV_per_block")){
    size_t EV_per_block=param_map.get_as_sizet("EV_per_block");
    if(EV_per_block<1){std::cerr<<"EV_per_block<1!"<<std::endl;exit(1);}
    nr_blocks=nr_print/EV_per_block;
    if((int)(nr_blocks*EV_per_block)<(int)nr_print)nr_blocks++;
  }

  std::cout<<"nr_blocks: "<<nr_blocks<<std::endl;
  Distribute_Block db_EV(nr_print,nr_blocks);


  for(size_t bl=0; bl<db_EV.get_nr_blocks(); bl++){
    size_t this_nr=db_EV.get_blocksize(bl);
    if(should_print){
      std::cout<<"block: "<<bl<<"/"<<db_EV.get_nr_blocks()<<std::endl;
    }

    // EV2[i] = H * EV[i] :
    std::vector<std::vector<std::complex<double> > > EV2(this_nr,  \
      std::vector<std::complex<double> >(hamil2.get_DIM(),0.));

    std::vector<double> eigenvalue(this_nr,0);
    std::vector<std::vector<std::complex<double> > >EV_local(this_nr, \
      std::vector<std::complex<double> >(db->get_local_blocksize(),0.));

    for(size_t i=0; i<this_nr; i++){  //For all eigenvectors:
      eigenvalue[i]=solver->get_eigenvalue(db_EV.get_start(bl)+i);
      solver->get_EV_local_part(db_EV.get_start(bl)+i, &EV_local[i][0]);
    }

    for(size_t j=0; j<db->get_local_blocksize(); j++){
      int j2=db->get_local_start()+j;
      std::vector< col_val > col_val_vec;
      Fermion_Double_Index fi=hamil->hilbert->get_Fermion_Double_Index(j2);
      hamil2.get_Matrix_row_block(col_val_vec, fi);
      for(size_t k=0; k<col_val_vec[0].size(); k++){
        int Hcol=col_val_vec[0][k].first;
        std::complex<double> Hval=std::conj(col_val_vec[0][k].second);
        for(size_t i=0; i<this_nr; i++){  
          EV2[i][Hcol]+=Hval*EV_local[i][j];
        }
      }
    }

    //Communicate EV2 to all threads
    if(should_print){
      std::cout<<"Communicating EV2..."<<std::endl;
    }
    for(size_t i=0; i<this_nr; i++){  
      Communicate_Complex_Vector_Template<Tsolver>(solver, EV2[i]);
    }

 
    //Calculate actual correction:
    if(should_print){
      std::cout<<"Calculating corrections..."<<std::endl;
    }

    std::vector<double> corr(this_nr,0.);
//    std::vector<double> norm_corr(this_nr,0.);
    std::vector<double> maxcontrib(this_nr,0.);
    std::vector<int> nr_small_denominator(this_nr,0);
 
   for(size_t j=db2->get_local_start(); j<db2->get_local_start()+db2->get_local_blocksize(); j++){
    //For those states that are not already in the smaller Hilbert space:
    Fermion_Double_Index fi=hamil2.hilbert->get_Fermion_Double_Index(j);
    int k=hamil->hilbert->get_index(fi);
    if(k<0){
      double e=hamil2.get_diagonal(fi);
      bool added=false;
      for(size_t i=0; i<this_nr; i++){  
        double numerator=(std::conj(EV2[i][j])*EV2[i][j]).real();
        double denominator=(e-eigenvalue[i]);

        if(print_new_states && fabs(denominator)<small_denominator_threshold){
          nr_small_denominator[i]++;
          newstates_indices.push_back(j);
          added=true;
          break;
        }
        corr[i]-=numerator/denominator;
//        norm_corr[i]+=numerator/(denominator*denominator);
        if(print_new_states){
          double contrib=abs(EV2[i][j]/denominator);
          if(contrib>maxcontrib[i])maxcontrib[i]=contrib;
          if(!added && contrib>threshold){
            newstates_indices.push_back(j);
            added=true;
          }
        }
      }
    }
   }
    
    //printing individual eigenvalues/corrections
    if(should_print){
      std::cout<<"Communicating and printing perturbative corrections..."<<std::endl;
    }
    for(size_t i=0; i<this_nr; i++){  
      //Communicate corr;
      {std::vector<std::complex<double> > c_corr(1,corr[i]); 
      Communicate_Complex_Vector_Template<Tsolver>(solver, c_corr);
      corr[i]=c_corr[0].real();}
      //Communicate norm_corr;
/*      {std::vector<std::complex<double> > c_norm_corr(1,norm_corr[i]); 
      Communicate_Complex_Vector_Template<Tsolver>(solver, c_norm_corr);
      norm_corr[i]=c_norm_corr[0].real();}
*/

      comm->Reduce_Max_Double(maxcontrib[i]);

      if(should_print){
        std::cout<<eigenvalue[i]<<" ";
        ofs<<eigenvalue[i]<<" ";
        if(nr_small_denominator[i]>0){
          std::cout<<"ILL_DEFINED_"<<nr_small_denominator[i]<<"_";
          ofs<<"ILL_DEFINED_"<<nr_small_denominator[i]<<"_";
        }
        std::cout<<corr[i];
//        std::cout<<" "<<norm_corr[i]
        std::cout<<std::endl;
        std::cout<<"Maximum contribution from new state: "<<maxcontrib[i]<<std::endl;
        ofs<<corr[i];
//        ofs<<" "<<norm_corr[i];
        ofs<<" "<<maxcontrib[i];
        ofs<<std::endl; 
      }
    }//<- loop over eigenvectors i
  }//<- loop over blocks


  //TODO: Communicate new states:
  if(should_print){
    std::cout<<"Communicating and printing new states..."<<std::endl;
  }
  Communicate_Append_Int_Vector_Template<Tsolver>(solver, newstates_indices);



  if(should_print)for(size_t j=0; j<newstates_indices.size(); j++){
    Fermion_Double_Index fi=hamil2.hilbert->get_Fermion_Double_Index(newstates_indices[j]);
    newstates.add_if_not_contained(fi);
  }


  //print:
  if(should_print){
    newstates.print(new_states_file);
    std::cout<<"Threshold: "<<threshold<<std::endl;
    std::cout<<"Adding "<<newstates.size()-existing_states_in_new_states_file;
    std::cout<<" states to new_states_file '"<<new_states_file;
    std::cout<<"' gives a total of "<<newstates.size()<<" states."<<std::endl;
    std::cout<<"Perturbative corrections done."<<std::endl;
  }

  delete comm;
  delete db2;
  delete db;
}



#endif