Commit ba112b4c authored by gbalduzz's avatar gbalduzz
Browse files

Fix combinatorial factor. Add parameter to decide if proposing partners on different sites or not.

parent c12e9245

include/dca/phys/dca_data/dca_data.hpp

+5 −0
Original line number Diff line number Diff line
@@ -232,6 +232,11 @@ public: // Optional members getters. 
              "non_density_interaction"));

    return *non_density_interactions_;

  }

  const auto& get_non_density_interactions() const {

    assert(non_density_interactions_);

    return *non_density_interactions_;

  }



  bool has_non_density_interactions() const {

    return (bool)non_density_interactions_;

  }

include/dca/phys/dca_step/cluster_solver/ctint/structs/interaction_vertices.hpp

+9 −5
Original line number Diff line number Diff line
@@ -40,6 +40,8 @@ struct InteractionElement { 
// Store the interaction terms and allow drawing a random vertex with strength |w|.

class InteractionVertices {

public:

  void initialize(double double_insertion_prob, bool all_sites_partnership);



  // Initialize vertices with a density-density Hamiltonian.

  // Precondition: Domain has the shape of dmn_variadic<Nu, Nu, Rdmn>

  template <class Nu, class Rdmn>
@@ -77,18 +79,20 @@ public: 
  }



  // Returns the number of possible partners for each non density-density interaction.

  int possiblePartners() const {

    const int partners = elements_.back().partners_id.size();

    // TODO: generalize if number of possible pairings is not constant or at the back.

    return partners;

  int possiblePartners(unsigned idx) const {

    assert(idx < elements_.size());

    return elements_[idx].partners_id.size();

  }



  std::vector<InteractionElement> elements_;



private:

  enum PartnershipType { NONE, SAME_SITE, ALL_SITES };



  std::vector<double> cumulative_weigths_;

  double total_weigth_ = 0;

  bool interband_propagator_ = false;

  PartnershipType partnership_type_ = NONE;

};



template <class Rng>
@@ -205,7 +209,7 @@ void InteractionVertices::checkForInterbandPropagators( 
  for (int r = 0; r < RDmn::dmn_size(); ++r)

    for (int b1 = 0; b1 < nb; ++b1)

      for (int b2 = 0; b2 < nb; ++b2) {

        if (r != r0 && b1 != b2 && std::abs(G_r_t(b1, 0, b2, 0, r, t0)) > 1e-8) {

        if (r != r0 && b1 != b2 && std::abs(G_r_t(b1, 0, b2, 0, r, t0)) > 1e-5) {

          interband_propagator_ = true;

          return;

        }

include/dca/phys/dca_step/cluster_solver/ctint/structs/solver_configuration.hpp

+3 −2
Original line number Diff line number Diff line
@@ -108,8 +108,9 @@ public: 
  // Return index corresponding to tag.

  int findTag(std::uint64_t tag) const;



  auto possiblePartners() const {

    return H_int_->possiblePartners();

  auto possiblePartners(unsigned idx) const {

    assert(idx < vertices_.size());

    return H_int_->possiblePartners(vertices_[idx].interaction_id);

  }



  friend io::Buffer& operator<<(io::Buffer& buff, const SolverConfiguration& config);

include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp

+4 −5
Original line number Diff line number Diff line
@@ -48,7 +48,6 @@ namespace ctint { 
template <linalg::DeviceType device_type, class Parameters, typename Real>

class CtintWalker;





template <class Parameters, typename Real = double>

class CtintWalkerBase {

public:
@@ -140,7 +139,7 @@ public: 


  static void setDMatrixAlpha(const std::array<double, 3>& alphas, bool adjust_dd);



  static void setInteractionVertices(Data& data);

  static void setInteractionVertices(const Data& data, const Parameters& parameters);



  float stealFLOPs() {

    auto flop = flop_;
@@ -178,7 +177,6 @@ protected: // Members. 


  const Real beta_;

  static inline constexpr int n_bands_ = Parameters::bands;

  const int possible_partners_;



  const Real total_interaction_;  // Space integrated interaction Hamiltonian.
@@ -220,7 +218,6 @@ CtintWalkerBase<Parameters, Real>::CtintWalkerBase(const Parameters& parameters_ 
                     parameters_.getDoubleUpdateProbability()),



      beta_(parameters_.get_beta()),

      possible_partners_(configuration_.possiblePartners()),

      total_interaction_(vertices_.integratedInteraction()) {}



template <class Parameters, typename Real>
@@ -338,8 +335,10 @@ void CtintWalkerBase<Parameters, Real>::setDMatrixAlpha(const std::array<double, 
}



template <class Parameters, typename Real>

void CtintWalkerBase<Parameters, Real>::setInteractionVertices(Data& data) {

void CtintWalkerBase<Parameters, Real>::setInteractionVertices(const Data& data,

                                                               const Parameters& parameters) {

  vertices_.reset();

  vertices_.initialize(parameters.getDoubleUpdateProbability(), parameters.getAllSitesPartnership());

  vertices_.initializeFromHamiltonian(data.H_interactions);

  if (data.has_non_density_interactions()) {

    vertices_.checkForInterbandPropagators(data.G0_r_t_cluster_excluded);

include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_cpu.hpp

+9 −8
Original line number Diff line number Diff line
@@ -86,7 +86,6 @@ protected: 
  using BaseClass::d_builder_ptr_;

  using BaseClass::total_interaction_;

  using BaseClass::beta_;

  using BaseClass::possible_partners_;

  using BaseClass::sign_;

  using BaseClass::thread_id_;
@@ -222,10 +221,11 @@ Real CtintWalker<linalg::CPU, Parameters, Real>::insertionProbability(const int 
  const Real combinatorial_factor =

      delta_vertices == 1 ? old_size + 1 : (old_size + 2) * (configuration_.nPartners(old_size) + 1);



  const Real strength_factor =

      delta_vertices == 1 ? -beta_ * total_interaction_ * configuration_.getSign(old_size)

  const Real strength_factor = delta_vertices == 1

                                   ? -beta_ * total_interaction_ * configuration_.getSign(old_size)

                                   : total_interaction_ * beta_ * beta_ *

                                std::abs(configuration_.getStrength(old_size)) * possible_partners_;

                                         std::abs(configuration_.getStrength(old_size)) *

                                         configuration_.possiblePartners(old_size);



  const Real det_ratio = det_ratio_[0] * det_ratio_[1];
@@ -265,8 +265,9 @@ Real CtintWalker<linalg::CPU, Parameters, Real>::removalProbability() { 
  const Real combinatorial_factor =

      n_removed == 1 ? n : n * configuration_.nPartners(removal_list_[0]);

  const Real strength_factor =

      n_removed == 1 ? -beta_ * total_interaction_ * configuration_.getSign(removal_list_[0])

                     : total_interaction_ * beta_ * beta_ * possible_partners_ *

      n_removed == 1

          ? -beta_ * total_interaction_ * configuration_.getSign(removal_list_[0])

          : total_interaction_ * beta_ * beta_ * configuration_.possiblePartners(removal_list_[0]) *

                std::abs(configuration_.getStrength(removal_list_[0]));



  const Real det_ratio = det_ratio_[0] * det_ratio_[1];