Commit 1d466648 authored by gbalduzz's avatar gbalduzz
Browse files

Store one configuration sample between dca runs, if store-configuration is true.

parent aa06c66a
Loading
Loading
Loading
Loading
+8 −0
Original line number Diff line number Diff line
@@ -62,6 +62,8 @@ public:

  std::string get_path();

  void erase(const std::string& name);

  template <typename arbitrary_struct_t>
  static void to_file(const arbitrary_struct_t& arbitrary_struct, const std::string& file_name);

@@ -114,6 +116,12 @@ public:
    return execute(name, static_cast<io::Buffer::Container>(buffer));
  }

  template <class T>
  void rewrite(const std::string& name, const T& obj) {
    erase(name);
    execute(name, obj);
  }

  operator bool() const {
    return static_cast<bool>(file_);
  }
+6 −6
Original line number Diff line number Diff line
@@ -28,7 +28,7 @@ class MPIPacking : public virtual MPIProcessorGrouping {
public:
  MPIPacking() {}

  template <typename scalar_type>
  template <typename scalar_type, typename = std::enable_if_t<std::is_scalar_v<scalar_type>>>
  int get_buffer_size(scalar_type item) const;
  template <typename scalar_type>
  int get_buffer_size(const std::basic_string<scalar_type>& str) const;
@@ -41,7 +41,7 @@ public:
  template <typename scalar_type, class dmn_type>
  int get_buffer_size(const func::function<scalar_type, dmn_type>& f) const;

  template <typename scalar_type>
  template <typename scalar_type, typename = std::enable_if_t<std::is_scalar_v<scalar_type>>>
  void pack(char* buffer, int size, int& off_set, scalar_type item) const;
  template <typename scalar_type>
  void pack(char* buffer, int size, int& off_set, const std::basic_string<scalar_type>& str) const;
@@ -57,7 +57,7 @@ public:
  void pack(char* buffer, int size, int& off_set,
            const func::function<scalar_type, dmn_type>& f) const;

  template <typename scalar_type>
  template <typename scalar_type, typename = std::enable_if_t<std::is_scalar_v<scalar_type>>>
  void unpack(char* buffer, int size, int& off_set, scalar_type& item) const;
  template <typename scalar_type>
  void unpack(char* buffer, int size, int& off_set, std::basic_string<scalar_type>& str) const;
@@ -72,7 +72,7 @@ public:
  void unpack(char* buffer, int size, int& off_set, func::function<scalar_type, dmn_type>& f) const;
};

template <typename scalar_type>
template <typename scalar_type, typename>
int MPIPacking::get_buffer_size(scalar_type /*item*/) const {
  int size(0);

@@ -175,7 +175,7 @@ int MPIPacking::get_buffer_size(const func::function<scalar_type, dmn_type>& f)
  return result;
}

template <typename scalar_type>
template <typename scalar_type, typename>
void MPIPacking::pack(char* buffer, int size, int& off_set, scalar_type item) const {
  const scalar_type* tPtr(&item);

@@ -263,7 +263,7 @@ void MPIPacking::pack(char* buffer, int size, int& off_set,
           MPIProcessorGrouping::get());
}

template <typename scalar_type>
template <typename scalar_type, typename>
void MPIPacking::unpack(char* buffer, int size, int& off_set, scalar_type& item) const {
  scalar_type tmp;

+0 −23
Original line number Diff line number Diff line
@@ -114,7 +114,6 @@ public:
  void initialize();
  void initializeH0_and_H_i();
  void initialize_G0();
  void initializeSigma(const std::string& filename);

  void compute_single_particle_properties();
  void compute_Sigma_bands();
@@ -529,28 +528,6 @@ void DcaData<Parameters>::initialize_G0() {
  G0_r_t_cluster_excluded = G0_r_t;
}

template <class Parameters>
void DcaData<Parameters>::initializeSigma(const std::string& filename) {
  if (concurrency_.id() == concurrency_.first()) {
    io::HDF5Reader reader;
    reader.open_file(filename);

    if (parameters_.adjust_chemical_potential()) {
      reader.open_group("parameters");
      reader.open_group("physics");
      reader.execute("chemical-potential", parameters_.get_chemical_potential());
      reader.close_group();
      reader.close_group();
    }

    reader.open_group("functions");
    reader.execute(Sigma);
    reader.close_group();
  }

  concurrency_.broadcast(parameters_.get_chemical_potential());
  concurrency_.broadcast(Sigma);
}

template <class Parameters>
void DcaData<Parameters>::compute_single_particle_properties() {
+109 −75
Original line number Diff line number Diff line
@@ -67,7 +67,7 @@ public:
  using HTS_solver_type =
      solver::HighTemperatureSeriesExpansionSolver<dca::linalg::CPU, ParametersType, DcaDataType>;

  DcaLoop(ParametersType& parameters_ref, DcaDataType& MOMS_ref, concurrency_type& concurrency_ref);
  DcaLoop(ParametersType& parameters_ref, DcaDataType& data__ref, concurrency_type& concurrency_ref);

  void write();

@@ -75,7 +75,7 @@ public:
  void execute();
  void finalize();

protected:
private:
  void adjust_chemical_potential();

  void perform_cluster_mapping();
@@ -95,11 +95,12 @@ protected:

  void logSelfEnergy(int i);

  ParametersType& parameters;
  DcaDataType& MOMS;
  concurrency_type& concurrency;
  void readInitialStatus(const std::string& filename);

  ParametersType& parameters_;
  DcaDataType& data_;
  concurrency_type& concurrency_;

private:
  DcaLoopData<ParametersType> DCA_info_struct;

  cluster_exclusion_type cluster_exclusion_obj;
@@ -121,28 +122,28 @@ protected:

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::DcaLoop(ParametersType& parameters_ref,
                                                                DcaDataType& MOMS_ref,
                                                                DcaDataType& data_ref,
                                                                concurrency_type& concurrency_ref)
    : parameters(parameters_ref),
      MOMS(MOMS_ref),
      concurrency(concurrency_ref),
    : parameters_(parameters_ref),
      data_(data_ref),
      concurrency_(concurrency_ref),

      DCA_info_struct(),

      cluster_exclusion_obj(parameters, MOMS),
      double_counting_correction_obj(parameters, MOMS),
      cluster_exclusion_obj(parameters_, data_),
      double_counting_correction_obj(parameters_, data_),

      cluster_mapping_obj(parameters),
      lattice_mapping_obj(parameters),
      cluster_mapping_obj(parameters_),
      lattice_mapping_obj(parameters_),

      update_chemical_potential_obj(parameters, MOMS, cluster_mapping_obj),
      update_chemical_potential_obj(parameters_, data_, cluster_mapping_obj),

      output_file_(concurrency.id() == concurrency.first() ? std::make_shared<io::HDF5Writer>(false)
      output_file_(concurrency_.id() == concurrency_.first() ? std::make_shared<io::HDF5Writer>(false)
                                                             : nullptr),

      monte_carlo_integrator_(parameters_ref, MOMS_ref, output_file_) {
  if (concurrency.id() == concurrency.first()) {
    file_name_ = parameters.get_directory() + parameters.get_filename_dca();
      monte_carlo_integrator_(parameters_ref, data_, output_file_) {
  if (concurrency_.id() == concurrency_.first()) {
    file_name_ = parameters_.get_directory() + parameters_.get_filename_dca();

    dca::util::SignalHandler::registerFile(output_file_);

@@ -152,13 +153,13 @@ DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::DcaLoop(ParametersType&

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::write() {
  if (concurrency.id() == concurrency.first()) {
  if (concurrency_.id() == concurrency_.first()) {
    std::cout << "\n\n\t\t start writing " << file_name_ << "\t" << dca::util::print_time() << "\n\n";

    output_file_->set_verbose(true);

    parameters.write(*output_file_);
    MOMS.write(*output_file_);
    parameters_.write(*output_file_);
    data_.write(*output_file_);
    monte_carlo_integrator_.write(*output_file_);
    DCA_info_struct.write(*output_file_);

@@ -177,31 +178,64 @@ template <typename ParametersType, typename DcaDataType, typename MCIntegratorTy
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::initialize() {
  int last_completed = -1;

  if (parameters.autoresume()) {  // Try to read state of previous run.
    last_completed = DCA_info_struct.tryToRead(file_name_ + ".tmp", concurrency);
  if (parameters_.autoresume()) {  // Try to read state of previous run.
    last_completed = DCA_info_struct.tryToRead(file_name_ + ".tmp", concurrency_);
  }
  if (last_completed >= 0) {
    if (concurrency.id() == concurrency.first())
    if (concurrency_.id() == concurrency_.first())
      std::cout << "\n   *******  Resuming DCA from iteration " << last_completed + 1 << "  *******\n"
                << std::endl;

    dca_iteration_ = std::min(last_completed + 1, parameters.get_dca_iterations() - 1);
    MOMS.initializeSigma(file_name_ + ".tmp");
    perform_lattice_mapping();
    dca_iteration_ = std::min(last_completed + 1, parameters_.get_dca_iterations() - 1);
    readInitialStatus(file_name_ + ".tmp");
  }
  else if (parameters.get_initial_self_energy() != "zero") {
    MOMS.initializeSigma(parameters.get_initial_self_energy());
    perform_lattice_mapping();
  else if (parameters_.get_initial_self_energy() != "zero") {
    readInitialStatus(parameters_.get_initial_self_energy());
  }

  if (concurrency_.id() == concurrency_.first()) {
    output_file_->open_file(file_name_ + ".tmp", parameters_.autoresume() ? false : true);
  }
}

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::readInitialStatus(
    const std::string& filename) {
  io::Buffer buffer;

  if (concurrency_.id() == concurrency_.first()) {
    io::HDF5Reader reader;
    reader.open_file(filename);

    if (parameters_.adjust_chemical_potential()) {
      reader.open_group("parameters");
      reader.open_group("physics");
      reader.execute("chemical-potential", parameters_.get_chemical_potential());
      reader.close_group();
      reader.close_group();
    }

  if (concurrency.id() == concurrency.first()) {
    output_file_->open_file(file_name_ + ".tmp", parameters.autoresume() ? false : true);
    reader.open_group("functions");
    reader.execute(data_.Sigma);
    reader.close_group();

    reader.open_group("Configurations");
    reader.execute("sample", buffer);
    reader.close_group();
  }

  concurrency_.broadcast(parameters_.get_chemical_potential());
  concurrency_.broadcast(data_.Sigma);
  concurrency_.broadcast(buffer);

  monte_carlo_integrator_.setSampleConfiguration(buffer);

  perform_lattice_mapping();
}

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::execute() {
  for (; dca_iteration_ < parameters.get_dca_iterations(); dca_iteration_++) {
  for (; dca_iteration_ < parameters_.get_dca_iterations(); dca_iteration_++) {
    adjust_chemical_potential();

    perform_cluster_mapping();
@@ -221,7 +255,7 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::execute() {
    logSelfEnergy(dca_iteration_);  // Write a check point.

    if (L2_Sigma_difference <
        parameters.get_dca_accuracy())  // set the acquired accuracy on |Sigma_QMC - Sigma_cg|
        parameters_.get_dca_accuracy())  // set the acquired accuracy on |Sigma_QMC - Sigma_cg|
      break;
  }
}
@@ -229,13 +263,13 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::execute() {
template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::finalize() {
  perform_cluster_mapping_self_energy();
  MOMS.compute_Sigma_bands();
  MOMS.compute_single_particle_properties();
  data_.compute_Sigma_bands();
  data_.compute_single_particle_properties();
}

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::adjust_chemical_potential() {
  if (parameters.adjust_chemical_potential())
  if (parameters_.adjust_chemical_potential())
    update_chemical_potential_obj.execute();
}

@@ -248,40 +282,40 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::perform_cluster_map

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::perform_cluster_mapping_self_energy() {
  if (concurrency.id() == concurrency.first())
  if (concurrency_.id() == concurrency_.first())
    std::cout << "\n\t\t coarsegrain-Selfenergy " << dca::util::print_time();

  profiler_type profiler("coarsegrain-Selfenergy", "DCA", __LINE__);

  if (parameters.do_dca_plus())
    cluster_mapping_obj.compute_S_K_w(MOMS.Sigma_lattice, MOMS.Sigma_cluster);
  if (parameters_.do_dca_plus())
    cluster_mapping_obj.compute_S_K_w(data_.Sigma_lattice, data_.Sigma_cluster);
  else
    MOMS.Sigma_cluster = MOMS.Sigma;
    data_.Sigma_cluster = data_.Sigma;

  MOMS.print_Sigma_QMC_versus_Sigma_cg();
  data_.print_Sigma_QMC_versus_Sigma_cg();

  symmetrize::execute<Lattice>(MOMS.Sigma_cluster, MOMS.H_symmetry);
  symmetrize::execute<Lattice>(data_.Sigma_cluster, data_.H_symmetry);
}

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::perform_cluster_mapping_Greens_function() {
  if (concurrency.id() == concurrency.first())
  if (concurrency_.id() == concurrency_.first())
    std::cout << "\n\t\t coarsegrain-Greens-function " << dca::util::print_time();

  profiler_type profiler("coarsegrain-Greens-function", "DCA", __LINE__);

  // Finite-size QMC
  if (parameters.do_finite_size_qmc())
    compute_G_k_w(MOMS.H_DCA, MOMS.Sigma, parameters.get_chemical_potential(),
                  parameters.get_coarsegraining_threads(), MOMS.G_k_w);
  if (parameters_.do_finite_size_qmc())
    compute_G_k_w(data_.H_DCA, data_.Sigma, parameters_.get_chemical_potential(),
                  parameters_.get_coarsegraining_threads(), data_.G_k_w);
  // DCA+
  else if (parameters.do_dca_plus())
    cluster_mapping_obj.compute_G_K_w(MOMS.Sigma_lattice, MOMS.G_k_w);
  else if (parameters_.do_dca_plus())
    cluster_mapping_obj.compute_G_K_w(data_.Sigma_lattice, data_.G_k_w);
  // Standard DCA
  else
    cluster_mapping_obj.compute_G_K_w(MOMS.Sigma, MOMS.G_k_w);
    cluster_mapping_obj.compute_G_K_w(data_.Sigma, data_.G_k_w);

  symmetrize::execute<Lattice>(MOMS.G_k_w, MOMS.H_symmetry);
  symmetrize::execute<Lattice>(data_.G_k_w, data_.H_symmetry);
}

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
@@ -296,7 +330,7 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::adjust_impurity_sel

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::perform_cluster_exclusion_step() {
  if (concurrency.id() == concurrency.first())
  if (concurrency_.id() == concurrency_.first())
    std::cout << "\n\t\t cluster-exclusion-step " << dca::util::print_time();

  profiler_type profiler("cluster-exclusion-step", "DCA", __LINE__);
@@ -328,25 +362,25 @@ template <typename ParametersType, typename DcaDataType, typename MCIntegratorTy
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::perform_lattice_mapping() {
  profiler_type profiler("lattice-mapping", "DCA", __LINE__);

  if (concurrency.id() == concurrency.first())
  if (concurrency_.id() == concurrency_.first())
    std::cout << "\n\t\t lattice-mapping " << dca::util::print_time();

  if (parameters.do_dca_plus() || parameters.doPostInterpolation()) {
    if (parameters.hts_approximation()) {
      DcaDataType MOMS_HTS(parameters);
  if (parameters_.do_dca_plus() || parameters_.doPostInterpolation()) {
    if (parameters_.hts_approximation()) {
      DcaDataType data__HTS(parameters_);

      MOMS_HTS.H_HOST = MOMS.H_HOST;
      MOMS_HTS.H_interactions = MOMS.H_interactions;
      data__HTS.H_HOST = data_.H_HOST;
      data__HTS.H_interactions = data_.H_interactions;

      HTS_solver_type HTS_solver(parameters, MOMS_HTS);
      HTS_solver_type HTS_solver(parameters_, data__HTS);

      lattice_mapping_obj.execute_with_HTS_approximation(
          MOMS_HTS, HTS_solver, cluster_mapping_obj, MOMS.Sigma, MOMS.Sigma_lattice_interpolated,
          MOMS.Sigma_lattice_coarsegrained, MOMS.Sigma_lattice);
          data__HTS, HTS_solver, cluster_mapping_obj, data_.Sigma, data_.Sigma_lattice_interpolated,
          data_.Sigma_lattice_coarsegrained, data_.Sigma_lattice);
    }
    else {
      lattice_mapping_obj.execute(MOMS.Sigma, MOMS.Sigma_lattice_interpolated,
                                  MOMS.Sigma_lattice_coarsegrained, MOMS.Sigma_lattice);
      lattice_mapping_obj.execute(data_.Sigma, data_.Sigma_lattice_interpolated,
                                  data_.Sigma_lattice_coarsegrained, data_.Sigma_lattice);
    }
  }
}
@@ -354,26 +388,26 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::perform_lattice_map
template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::update_DCA_loop_data_functions(int i) {
  DCA_info_struct.density(i) = update_chemical_potential_obj.compute_density();
  DCA_info_struct.chemical_potential(i) = parameters.get_chemical_potential();
  DCA_info_struct.chemical_potential(i) = parameters_.get_chemical_potential();

  if (concurrency.id() == concurrency.first()) {
  if (concurrency_.id() == concurrency_.first()) {
    std::cout << "\n\n\t\t\t total-density : " << DCA_info_struct.density(i)
              << "\t (time : " << dca::util::print_time() << ")\n\n";
  }

  for (int l1 = 0; l1 < b::dmn_size() * s::dmn_size(); l1++)
    DCA_info_struct.orbital_occupancies(l1, i) = MOMS.orbital_occupancy(l1);
    DCA_info_struct.orbital_occupancies(l1, i) = data_.orbital_occupancy(l1);

  for (int l1 = 0; l1 < b::dmn_size() * s::dmn_size(); l1++)
    for (int k_ind = 0; k_ind < k_DCA::dmn_size(); k_ind++)
      DCA_info_struct.n_k(l1, k_ind, i) = 1. - std::abs(MOMS.G_k_t(l1, l1, k_ind, 0));
      DCA_info_struct.n_k(l1, k_ind, i) = 1. - std::abs(data_.G_k_t(l1, l1, k_ind, 0));

  for (int l1 = 0; l1 < b::dmn_size() * s::dmn_size(); l1++)
    for (int k_ind = 0; k_ind < k_DCA::dmn_size(); k_ind++)
      // TODO: Use t::dmn_size() instead of parameters.get_sp_time_intervals().
      // TODO: Use t::dmn_size() instead of parameters_.get_sp_time_intervals().
      DCA_info_struct.A_k(l1, k_ind, i) =
          std::abs(MOMS.G_k_t(l1, l1, k_ind, parameters.get_sp_time_intervals() / 2)) *
          parameters.get_beta() / M_PI;
          std::abs(data_.G_k_t(l1, l1, k_ind, parameters_.get_sp_time_intervals() / 2)) *
          parameters_.get_beta() / M_PI;
}

template <typename ParametersType, typename DcaDataType, typename MCIntegratorType>
@@ -382,12 +416,12 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::logSelfEnergy(int i

  if (output_file_) {
    output_file_->open_group("functions");
    output_file_->execute(MOMS.Sigma);
    output_file_->execute(data_.Sigma);
    output_file_->close_group();

    output_file_->open_group("parameters");
    output_file_->open_group("physics");
    output_file_->execute("chemical-potential", parameters.get_chemical_potential());
    output_file_->execute("chemical-potential", parameters_.get_chemical_potential());
    output_file_->close_group();
    output_file_->close_group();

@@ -395,7 +429,7 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType>::logSelfEnergy(int i
  }
}

}  // namespace phys
}  // namespace dca
}  // namespace dca

#endif  // DCA_PHYS_DCA_LOOP_DCA_LOOP_HPP
+2 −0
Original line number Diff line number Diff line
@@ -100,6 +100,8 @@ public:
  // Precondition: The accumulator_ data has not been averaged, i.e. finalize has not been called.
  auto local_G_k_w() const;

  void setSampleConfiguration(const io::Buffer&) {}

protected:
  void warmUp(Walker& walker);

Loading