Unverified Commit cce8d0b2 authored by Peter Doak's avatar Peter Doak Committed by GitHub
Browse files

Merge pull request #143 from gbalduzz/fix_develop_137

Properly store multiple G4 channel in the parameters.
parents 0067de9b 15c346d6
Loading
Loading
Loading
Loading
+20 −13
Original line number Diff line number Diff line
@@ -50,17 +50,26 @@ public:
  function(const std::string& name = default_name_);

  // Copy constructor
  // Constructs the function with the name name and a copy of the elements of other.
  // Constructs the function with the a copy of elements and name of other.
  // Precondition: The other function has been resetted, if the domain had been initialized after
  //               the other function's construction.
  function(const function<scalartype, domain>& other, const std::string& name = default_name_);
  function(const function<scalartype, domain>& other);
  // Same as above, but with name = 'name'.
  function(const function<scalartype, domain>& other, const std::string& name) : function(other) {
    name_ = name;
  }

  // Move constructor
  // Constructs the function with the name name and the elements of other using move semantics.
  // Constructs the function with elements and name of other using move semantics.
  // Precondition: The other function has been resetted, if the domain had been initialized after
  //               the other function's construction.
  // Postcondition: The other function is in a non-specified state.
  function(function<scalartype, domain>&& other, const std::string& name = default_name_);
  function(function<scalartype, domain>&& other);
  // Same as above, but with name = 'name'.
  function(function<scalartype, domain>&& other, const std::string& name)
      : function(std::move(other)) {
    name_ = name;
  }

  // Copy assignment operator
  // Replaces the function's elements with a copy of the elements of other.
@@ -191,7 +200,6 @@ public:
  void operator*=(const function<scalartype, domain>& other);
  void operator/=(const function<scalartype, domain>& other);


  void operator=(scalartype c);
  void operator+=(scalartype c);
  void operator-=(scalartype c);
@@ -267,9 +275,8 @@ function<scalartype, domain>::function(const std::string& name)
}

template <typename scalartype, class domain>
function<scalartype, domain>::function(const function<scalartype, domain>& other,
                                       const std::string& name)
    : name_(name),
function<scalartype, domain>::function(const function<scalartype, domain>& other)
    : name_(other.name_),
      function_type(__PRETTY_FUNCTION__),
      dmn(),
      Nb_elements(dmn.get_size()),
@@ -286,8 +293,8 @@ function<scalartype, domain>::function(const function<scalartype, domain>& other
}

template <typename scalartype, class domain>
function<scalartype, domain>::function(function<scalartype, domain>&& other, const std::string& name)
    : name_(name),
function<scalartype, domain>::function(function<scalartype, domain>&& other)
    : name_(std::move(other.name_)),
      function_type(__PRETTY_FUNCTION__),
      dmn(),
      Nb_elements(dmn.get_size()),
@@ -623,7 +630,7 @@ void function<scalartype, domain>::unpack(const concurrency_t& concurrency, char
  concurrency.unpack(buffer, buffer_size, position, *this);
}

}  // func
}  // dca
}  // namespace func
}  // namespace dca

#endif  // DCA_FUNCTION_FUNCTION_HPP
+2 −1
Original line number Diff line number Diff line
@@ -233,7 +233,8 @@ void BseClusterSolver<ParametersType, DcaDataType, ScalarType>::load_G_II_0(
      for (int n2 = 0; n2 < b::dmn_size(); n2++) {
        for (int m1 = 0; m1 < b::dmn_size(); m1++) {
          for (int m2 = 0; m2 < b::dmn_size(); m2++) {
            switch (parameters.get_four_point_type()) {
            // TODO: allow more than one channel.
            switch (parameters.get_four_point_channels()[0]) {
              case PARTICLE_HOLE_TRANSVERSE: {
                G_II_0(n1, n2, k, w_vertex, m1, m2, k, w_vertex) =
                    -data_.G_k_w(n1, e_UP, m2, e_UP, k, w) *
+6 −6
Original line number Diff line number Diff line
@@ -401,7 +401,8 @@ void BseLatticeSolver<ParametersType, DcaDataType, ScalarType>::diagonalizeGamma
#ifndef DCA_ANALYSIS_TEST_WITH_FULL_DIAGONALIZATION
    // Diagonalize the symmetric matrix \sqrt{\chi_0}\Gamma\sqrt{\chi_0}.
    // The origin in momentum space has always index = 0.
    if (parameters.get_four_point_type() == PARTICLE_PARTICLE_UP_DOWN &&
    // TODO: loop over multiple channels.
    if (parameters.get_four_point_channels()[0] == PARTICLE_PARTICLE_UP_DOWN &&
        parameters.get_four_point_momentum_transfer_index() == 0 &&
        parameters.get_four_point_frequency_transfer() == 0) {
      diagonalizeGammaChi0Symmetric();
@@ -618,8 +619,7 @@ void BseLatticeSolver<ParametersType, DcaDataType, ScalarType>::diagonalize_fold

  {
    if (concurrency.id() == concurrency.first())
      std::cout << "\n\n\t diagonalize P_Gamma_chi_0_lattice_P " << dca::util::print_time()
                << " ...";
      std::cout << "\n\n\t diagonalize P_Gamma_chi_0_lattice_P " << dca::util::print_time() << " ...";

    dca::linalg::Vector<std::complex<ScalarType>, dca::linalg::CPU> L("L (BseLatticeSolver)", M);

@@ -924,8 +924,8 @@ void BseLatticeSolver<ParametersType, DcaDataType, ScalarType>::characterizeLead
  }
}

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

#endif  // DCA_PHYS_DCA_ANALYSIS_BSE_SOLVER_BSE_LATTICE_SOLVER_HPP
+7 −3
Original line number Diff line number Diff line
@@ -101,6 +101,10 @@ BseSolver<ParametersType, DcaDataType>::BseSolver(ParametersType& parameters, Dc

      wave_functions_names_("wave-functions-names"),
      harmonics_("harmonics") {
  if (parameters_.get_four_point_channels().size() != 1) {
    throw(std::logic_error("The analysis code can analyze exactly one channel."));
  }

  initialize_wave_functions();

  {
@@ -227,8 +231,8 @@ void BseSolver<ParametersType, DcaDataType>::calculateSusceptibilities() {
  bse_lattice_solver_.diagonalizeGammaChi0();
}

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

#endif  // DCA_PHYS_DCA_ANALYSIS_BSE_SOLVER_BSE_SOLVER_HPP
+13 −57
Original line number Diff line number Diff line
@@ -287,60 +287,11 @@ DcaData<Parameters>::DcaData(Parameters& parameters_ref)
  // Reserve storage in advance such that we don't have to copy elements when we fill the vector.
  // We want to avoid copies because function's copy ctor does not copy the name (and because copies
  // are expensive).
  G4_.reserve(parameters_.numG4Channels());

  for (auto channel : parameters_.get_four_point_channels()) {
    // Allocate memory for G4.
  // Ensure backward compatibility.
  if (parameters_.get_four_point_type() != NONE)
    G4_.emplace_back("G4");

  // Check which four point types to accumulate.
  else {
    if (parameters_.accumulateG4ParticleHoleTransverse())
      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_TRANSVERSE));

    if (parameters_.accumulateG4ParticleHoleMagnetic())
      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_MAGNETIC));

    if (parameters_.accumulateG4ParticleHoleCharge())
      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_CHARGE));

    if (parameters_.accumulateG4ParticleHoleLongitudinalUpUp())
      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_UP));

    if (parameters_.accumulateG4ParticleHoleLongitudinalUpDown())
      G4_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_DOWN));

    if (parameters_.accumulateG4ParticleParticleUpDown())
      G4_.emplace_back("G4_" + toString(PARTICLE_PARTICLE_UP_DOWN));
  }

    G4_.emplace_back("G4_" + toString(channel));
    // Allocate memory for error on G4.
  if (parameters_.get_error_computation_type() != ErrorComputationType::NONE) {
    G4_err_.reserve(parameters_.numG4Channels());

    if (parameters_.get_four_point_type() != NONE)
      G4_err_.emplace_back("G4-error");

    else {
      if (parameters_.accumulateG4ParticleHoleTransverse())
        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_TRANSVERSE) + "_err");

      if (parameters_.accumulateG4ParticleHoleMagnetic())
        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_MAGNETIC) + "_err");

      if (parameters_.accumulateG4ParticleHoleCharge())
        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_CHARGE) + "_err");

      if (parameters_.accumulateG4ParticleHoleLongitudinalUpUp())
        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_UP) + "_err");

      if (parameters_.accumulateG4ParticleHoleLongitudinalUpDown())
        G4_err_.emplace_back("G4_" + toString(PARTICLE_HOLE_LONGITUDINAL_UP_DOWN) + "_err");

      if (parameters_.accumulateG4ParticleParticleUpDown())
        G4_err_.emplace_back("G4_" + toString(PARTICLE_PARTICLE_UP_DOWN) + "_err");
    }
    G4_err_.emplace_back("G4_" + toString(channel) + "_err");
  }
}

@@ -374,7 +325,7 @@ void DcaData<Parameters>::read(std::string filename) {

  concurrency_.broadcast_object(Sigma);

  if (parameters_.accumulateG4()) {
  if (parameters_.isAccumulatingG4()) {
    concurrency_.broadcast_object(G_k_w);

    for (auto& G4_channel : G4_)
@@ -397,9 +348,14 @@ void DcaData<Parameters>::read(Reader& reader) {

  reader.execute(Sigma);

  if (parameters_.accumulateG4()) {
  if (parameters_.isAccumulatingG4()) {
    reader.execute(G_k_w);

    // Try to read G4 with a legacy name.
    if (parameters_.get_four_point_channels().size() == 1) {
      reader.execute("G4", G4_[0]);
    }

    for (auto& G4_channel : G4_)
      reader.execute(G4_channel);
  }
@@ -498,7 +454,7 @@ void DcaData<Parameters>::write(Writer& writer) {
    writer.execute(G0_r_t_cluster_excluded);
  }

  if (parameters_.accumulateG4()) {
  if (parameters_.isAccumulatingG4()) {
    if (!(parameters_.dump_cluster_Greens_functions())) {
      writer.execute(G_k_w);
      writer.execute(G_k_w_err_);
Loading