Commit 2587c456 authored by gbalduzz's avatar gbalduzz
Browse files

Better handling of double counting.

parent 8f266b7e

include/dca/phys/dca_data/dca_data.hpp

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

    return *non_density_interactions_;

  }

  const auto& get_nondensity_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/ctint_cluster_solver.hpp

+1 −1
Original line number Diff line number Diff line
@@ -153,7 +153,7 @@ CtintClusterSolver<device_t, Parameters, use_submatrix>::CtintClusterSolver( 
  Walker::setDMatrixBuilder(g0_, RDmn::parameter_type::get_subtract_matrix(),

                            label_dmn_.get_branch_domain_steps(), parameters_.getAlphas());



  Walker::setInteractionVertices(parameters_, data_);

    Walker::setInteractionVertices(data_);



  if (concurrency_.id() == concurrency_.first())

    std::cout << "\n\n\t CT-INT Integrator is born \n\n";

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

+38 −10
Original line number Diff line number Diff line
@@ -42,8 +42,7 @@ public: 
  // Initialize vertices with a density-density Hamiltonian.

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

  template <class Nu, class Rdmn>

  void initializeFromHamiltonian(const func::function<double, func::dmn_variadic<Nu, Nu, Rdmn>>& H_int,

                                 bool double_counted = true);

  void initializeFromHamiltonian(const func::function<double, func::dmn_variadic<Nu, Nu, Rdmn>>& H_int);

  template <class Nu, class Rdmn>

  void initializeFromNonDensityHamiltonian(

      const func::function<double, func::dmn_variadic<Nu, Nu, Nu, Nu, Rdmn>>& H_int);
@@ -78,16 +77,45 @@ private: 


template <class Nu, class Rdmn>

void InteractionVertices::initializeFromHamiltonian(

    const func::function<double, func::dmn_variadic<Nu, Nu, Rdmn>>& H_int, bool double_counted) {

    const func::function<double, func::dmn_variadic<Nu, Nu, Rdmn>>& H_int) {

  // Note: Use this version of the code if double counting in the interaction hamiltonian is removed.

  //

  //  for (ushort nu1 = 0; nu1 < Nu::dmn_size(); nu1++) {

  //    for (ushort nu2 = 0; nu2 < Nu::dmn_size(); nu2++)

  //      for (ushort delta_r = 0; delta_r < Rdmn::dmn_size(); delta_r++) {

  //        const double value = H_int(nu1, nu2, delta_r);

  //        if (std::abs(value) < 1e-8)

  //          continue;

  //        for (ushort r1 = 0; r1 < Rdmn::dmn_size(); r1++) {

  //          const ushort r2 = Rdmn::parameter_type::subtract(delta_r, r1); // delta_r = r1 - r2

  //          insertElement(InteractionElement{{r1, r1, r2, r2}, {nu1, nu1, nu2, nu2}, value, short(-1)});

  //        }

  //      }

  //  }



  // Assume the density-density interaction Hamiltonian function is double counted, i.e.

  // H(b1, b2, r1 - r2) == H(b2, b1, r -r2) and both terms describe the same addendum to the

  // physical Hamiltonian.

  func::function<bool, func::dmn_variadic<Nu, Nu, Rdmn>> already_inserted;

  const int r0 = Rdmn::parameter_type::origin_index();

  for (ushort nu1 = 0; nu1 < Nu::dmn_size(); nu1++) {

    const ushort nu2_start = double_counted ? nu1 + 1 : 0;

    for (ushort nu2 = nu2_start; nu2 < Nu::dmn_size(); nu2++)

    for (ushort nu2 = 0; nu2 < Nu::dmn_size(); nu2++)

      for (ushort delta_r = 0; delta_r < Rdmn::dmn_size(); delta_r++) {

        const double value = H_int(nu1, nu2, delta_r);

        if (std::abs(value) < 1e-8)

          continue;



        const double minus_delta_r = Rdmn::parameter_type::subtract(delta_r, r0);

        if (std::abs(H_int(nu1, nu2, delta_r) - H_int(nu2, nu1, minus_delta_r)) > 1e-8)

          throw(std::logic_error("The interaction double counting is not consistent."));



        if (already_inserted(nu2, nu1, minus_delta_r))  // Avoid double counting.

          continue;



        // Insert

        already_inserted(nu1, nu2, delta_r) = true;

        for (ushort r1 = 0; r1 < Rdmn::dmn_size(); r1++) {

          const ushort r2 = Rdmn::parameter_type::add(delta_r, r1);

          const ushort r2 = Rdmn::parameter_type::subtract(delta_r, r1);  // delta_r = r1 - r2

          insertElement(InteractionElement{{r1, r1, r2, r2}, {nu1, nu1, nu2, nu2}, value, short(-1)});

        }

      }
@@ -122,9 +150,9 @@ void InteractionVertices::initializeFromNonDensityHamiltonian( 
          }

}



}  // ctint

}  // solver

}  // phys

}  // dca

}  // namespace ctint

}  // namespace solver

}  // namespace phys

}  // namespace dca



#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_INTERACTION_VERTICES_HPP

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

+4 −4
Original line number Diff line number Diff line
@@ -142,7 +142,7 @@ public: 
                                const std::vector<std::size_t>& sbdm_step,

                                const std::array<double, 3>& alphas);



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

  static void setInteractionVertices(Data &data);



  float stealFLOPs() {

    auto flop = flop_;
@@ -326,9 +326,9 @@ void CtintWalkerBase<Parameters>::setDMatrixBuilder( 
}



template <class Parameters>

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

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

  vertices_.reset();

  vertices_.initializeFromHamiltonian(data.H_interactions, parameters.doubleCountedInteraction());

  vertices_.initializeFromHamiltonian(data.H_interactions);

  if (data.has_non_density_interactions())

    vertices_.initializeFromNonDensityHamiltonian(data.get_non_density_interactions());

}

include/dca/phys/models/traits.hpp

+3 −3
Original line number Diff line number Diff line
@@ -54,8 +54,8 @@ template <class Lattice, class HType, class Parameters> 
typename std::enable_if_t<not has_non_density_interaction<Lattice>::value, void> initializeNonDensityInteraction(

    HType& /*interaction*/, const Parameters& /*pars*/) {}



}  // models

}  // phys

}  // dca

}  // namespace models

}  // namespace phys

}  // namespace dca



#endif  // DCA_PHYS_MODELS_TRAITS_HPP