Candidate for ComplexG0 merge

CMakeLists.txt

@@ -169,7 +169,6 @@ set(DCA_LIBS
   quantum_domains
   time_and_frequency_domains
   signals
-  symmetrization
   coarsegraining
   ${DCA_CONCURRENCY_LIB}
   ${DCA_THREADING_LIBS}

applications/analysis/chi_q_omega.cpp

@@ -14,6 +14,7 @@
 #include <string>
 
 // Defines Concurrency, Threading, ParametersType, DcaData, and BseSolver.
+#include "dca/config/mc_options.hpp"
 #include "dca/config/analysis.hpp"
 #include "dca/config/cmake_options.hpp"
 #include "dca/io/json/json_reader.hpp"

applications/analysis/main_analysis.cpp

@@ -13,7 +13,8 @@
 
 #include <iostream>
 #include <string>
-
+// Defined Precision and Scalar Type
+#include "dca/config/mc_options.hpp"
 // Defines Concurrency, Threading, ParametersType, DcaData, and BseSolver.
 #include "dca/config/analysis.hpp"
 #include "dca/config/cmake_options.hpp"
@@ -63,7 +64,6 @@ int main(int argc, char** argv) {
   adios2::ADIOS adios;
 
   if (dca::io::stringToIOType(parameters.get_output_format()) == dca::io::IOType::ADIOS2) {
-    int rank = concurrency.id();
     std::cout << "\nProcessor " << concurrency.id() << " is writing data." << std::endl;
     dca::io::Writer writer(adios, concurrency, parameters.get_output_format(), true);
     std::string filename_bse(parameters.get_directory() + parameters.getAppropriateFilenameAnalysis());

applications/dca/CMakeLists.txt

@@ -5,7 +5,7 @@ if (DCA_BUILD_DCA)
   target_include_directories(main_dca PRIVATE ${DCA_INCLUDE_DIRS})
 
   if (DCA_HAVE_GPU)
-    target_link_libraries(main_dca PRIVATE ${DCA_GPU_LIBS})
+    target_link_libraries(main_dca PRIVATE ${DCA_GPU_LIBS} mc_kernels g0_interpolation mc_kernels ctint mc_kernels)
   endif()
 
   target_link_libraries(main_dca PUBLIC FFTW::Double signals  ${DCA_LIBS})

applications/dca/main_dca.cpp

@@ -64,8 +64,14 @@ int dca_main(int argc, char** argv) {
     // Create the parameters object from the input file.
     ParametersType parameters(dca::util::GitVersion::string(), concurrency);
     parameters.read_input_and_broadcast<dca::io::JSONReader>(input_file);
+    if(concurrency.id() == concurrency.first())
+      std::cout << "Input read and broadcast.\n";
     parameters.update_model();
+    if(concurrency.id() == concurrency.first())
+      std::cout << "Model updated.\n";
     parameters.update_domains();
+    if(concurrency.id() == concurrency.first())
+      std::cout << "Domains updated.\n";
 
     dca::DistType distribution = parameters.get_g4_distribution();
     {

build-aux/summit_load_modules.sh

@@ -6,9 +6,9 @@
 # Usage: source summit_load_modules.sh
 
 module reset
-module load gcc/11.2.0
-module load cuda/11.5.2 # ldd shows magma is built with this cuda
-module load magma/2.6.2
+module load gcc/10.2.0
+module load cuda/11.1.1 # ldd shows magma is built with this cuda
+module load magma/2.6.1
 module load hdf5
 module load fftw
 module load cmake/3.21.3

cmake/dca_config.cmake

@@ -121,9 +121,9 @@ endif()
 
 # Lattice type
 set(DCA_LATTICE "square" CACHE STRING "Lattice type, options are: bilayer | square | triangular |
-    Kagome | hund | twoband_Cu | threeband | FeAs | material_NiO | material_FeSn ")
+    Kagome | hund | twoband_Cu | threeband | Rashba_Hubbard | Moire_Hubbard | FeAs | material_NiO | material_FeSn ")
 set_property(CACHE DCA_LATTICE PROPERTY STRINGS bilayer square triangular Kagome hund twoband_Cu threeband
-             FeAs material_NiO material_FeSn)
+             Rashba_Hubbard Moire_Hubbard FeAs material_NiO material_FeSn)
 
 if (DCA_LATTICE STREQUAL "bilayer")
   set(DCA_LATTICE_TYPE dca::phys::models::bilayer_lattice<PointGroup>)
@@ -145,12 +145,18 @@ elseif (DCA_LATTICE STREQUAL "Kagome")
     "dca/phys/models/analytic_hamiltonians/Kagome_hubbard.hpp")
 elseif (DCA_LATTICE STREQUAL "hund")
   set(DCA_LATTICE_TYPE dca::phys::models::HundLattice<PointGroup>)
-
 elseif (DCA_LATTICE STREQUAL "threeband")
   set(DCA_LATTICE_TYPE dca::phys::models::ThreebandHubbard<PointGroup>)
   set(DCA_LATTICE_INCLUDE
     "dca/phys/models/analytic_hamiltonians/threeband_hubbard.hpp")
-
+elseif (DCA_LATTICE STREQUAL "Rashba_Hubbard")
+  set(DCA_LATTICE_TYPE dca::phys::models::RashbaHubbard<PointGroup>)
+  set(DCA_LATTICE_INCLUDE
+    "dca/phys/models/analytic_hamiltonians/rashba_hubbard.hpp")
+elseif (DCA_LATTICE STREQUAL "Moire_Hubbard")
+  set(DCA_LATTICE_TYPE dca::phys::models::moire_hubbard<PointGroup>)
+  set(DCA_LATTICE_INCLUDE
+    "dca/phys/models/analytic_hamiltonians/Moire_Hubbard.hpp")
 elseif (DCA_LATTICE STREQUAL "twoband_chain")
   set(DCA_LATTICE_TYPE dca::phys::models::twoband_chain<dca::phys::domains::no_symmetry<1>>)
   set(DCA_LATTICE_INCLUDE
@@ -174,7 +180,7 @@ elseif (DCA_LATTICE STREQUAL "material_FeSn")
       "dca/phys/models/material_hamiltonians/material_lattice.hpp")
   set(DCA_MODEL_IS_MATERIAL_LATTICE ON CACHE BOOL "is the model a material lattice")
 else()
-  message(FATAL_ERROR "Please set DCA_LATTICE to a valid option: bilayer | square | triangular | Kagome | hund | twoband_Cu | threeband | FeAs | material_NiO | material_FeSn.")
+  message(FATAL_ERROR "Please set DCA_LATTICE to a valid option: bilayer | square | triangular | Kagome | hund | twoband_Cu | threeband | Rashba_Hubbard | Moire_Hubbard | FeAs | material_NiO | material_FeSn.")
 endif()
 
 # Model type
@@ -390,15 +396,15 @@ option(DCA_WITH_SINGLE_PRECISION_TP_MEASUREMENTS "Measure two particle function
 
 if (DCA_WITH_SINGLE_PRECISION_MC)
   set(DCA_WITH_SINGLE_PRECISION_TP_MEASUREMENTS ON CACHE BOOL "Measure two particle function in single precision." FORCE)
-  set(MC_SCALAR float)
+  set(MC_REAL float)
 else()
-  set(MC_SCALAR double)
+  set(MC_REAL double)
 endif()
 
 if (DCA_WITH_SINGLE_PRECISION_TP_MEASUREMENTS)
-  set(TP_ACCUMULATION_SCALAR float)
+  set(TP_ACCUMULATION_PRECISION float)
 else()
-  set(TP_ACCUMULATION_SCALAR double)
+  set(TP_ACCUMULATION_PRECISION double)
 endif()
 
 

cmake/dca_testing.cmake

@@ -184,7 +184,7 @@ function(dca_add_gtest name)
 
     add_test(NAME ${name}
              COMMAND ${TEST_RUNNER} ${MPIEXEC_NUMPROC_FLAG} ${DCA_ADD_GTEST_MPI_NUMPROC}
-                     ${MPIEXEC_PREFLAGS} ${SMPIARGS_FLAG_MPI} ${CVD_LAUNCHER} "$<TARGET_FILE:${name}>"
+                     ${MPIEXEC_PREFLAGS} ${SMPIARGS_FLAG_MPI} ${MPIEXEC_POSTFLAGS} ${CVD_LAUNCHER} "$<TARGET_FILE:${name}>"
                      ${DCA_TESTING_FLAGS})
                  target_link_libraries(${name} PRIVATE ${MPI_C_LIBRARIES})
   else()

examples/NiO_CTINT/input_NiO_ctint.json

@@ -37,7 +37,7 @@
 
   "DCA" : {
     "accuracy" : 0,
-      "interacting-orbitals" : [ 0, 1, 2, 3, 4, 5, 6, 7 ],
+    "interacting-orbitals" : [ 0, 1, 2, 3, 4, 5, 6, 7 ],
     "iterations" : 9,
     "self-energy-mixing-factor" : 0.75,
 

include/dca/application/dca_loop_dispatch.hpp

@@ -25,12 +25,16 @@ class DCALoopDispatch {
     // Create and initialize the DCA data object.
     DcaDataType<DT> dca_data(parameters);
     dca_data.initialize();
+    if(concurrency.id() == concurrency.first())
+      std::cout << "dca_data initialized.\n";
     DcaLoopType<DT> dca_loop(parameters, dca_data, concurrency);
     {
       Profiler profiler(__FUNCTION__, __FILE__, __LINE__);
 
       try {
         dca_loop.initialize();
+	if(concurrency.id() == concurrency.first())
+	  std::cout << "dca_loop initialized.\n";
       }
       catch (const std::exception& exc) {
         std::cout << "unhandled exception in dca_loop.initialize(): " << exc.what() << std::endl;
@@ -56,9 +60,10 @@ class DCALoopDispatch {
       //      if (concurrency.id() == concurrency.first()) {
       // std::cout << "\nProcessor " << concurrency.id() << " is writing data." << std::endl;
       dca_loop.write();
-      std::cout << "final output written to file: " << parameters.get_filename_dca() << '\n';
-
-      std::cout << "\nFinish time: " << dca::util::print_time() << "\n" << std::endl;
+      if(concurrency.id() == concurrency.first()) {
+	std::cout << "final output written to file: " << parameters.get_filename_dca() << '\n';
+	std::cout << "\nFinish time: " << dca::util::print_time() << "\n" << std::endl;
+      }
       //      }
     }
   }

include/dca/config/analysis.hpp.in

@@ -24,8 +24,11 @@
 constexpr dca::ClusterSolverId solver_name = @DCA_CLUSTER_SOLVER_NAME@;
 
 #include "dca/phys/parameters/parameters.hpp"
-using ParametersType = dca::phys::params::Parameters<Concurrency, Threading, Profiler, Model,
-                                                     void /*RandomNumberGenerator*/, solver_name>;
+using ParametersType = dca::phys::params::Parameters<
+    Concurrency, Threading, Profiler, Model, void, solver_name,
+    dca::NumericalTraits<dca::config::McOptions::MC_REAL,
+                         typename dca::util::ScalarSelect<dca::config::McOptions::MC_REAL,
+                                                          Model::lattice_type::complex_g0>::type>>;
 
 #include "dca/phys/dca_data/dca_data.hpp"
 using DcaDataType = dca::phys::DcaData<ParametersType>;

include/dca/config/dca.hpp.in

@@ -25,32 +25,37 @@
 #include "dca/phys/dca_step/cluster_solver/cluster_solver_id.hpp"
 constexpr dca::ClusterSolverId solver_name = @DCA_CLUSTER_SOLVER_NAME@;
 
+#include "dca/config/mc_options.hpp"
+
 #include "dca/phys/parameters/parameters.hpp"
-using ParametersType = dca::phys::params::Parameters<Concurrency, Threading, Profiler, Model,
-                                                     RandomNumberGenerator, solver_name>;
+using ParametersType = dca::phys::params::Parameters<
+    Concurrency, Threading, Profiler, Model, RandomNumberGenerator, solver_name,
+    dca::NumericalTraits<dca::config::McOptions::MC_REAL,
+                         typename dca::util::ScalarSelect<dca::config::McOptions::MC_REAL,
+                                                          Model::lattice_type::complex_g0>::type>>;
 
 #include "dca/phys/dca_data/dca_data.hpp"
-template<dca::DistType DIST>
+template <dca::DistType DIST>
 using DcaDataType = dca::phys::DcaData<ParametersType, DIST>;
 
 // You know whats probably worse than using the C preprocessor,
 // generating code using cmake.
 #include "@DCA_CLUSTER_SOLVER_INCLUDE@"
-template<dca::DistType DIST>
+template <dca::DistType DIST>
 using ClusterSolverBaseType = @DCA_CLUSTER_SOLVER_TYPE@;
 
 // TODO: Clean this up (remove #ifdef).
 #ifdef DCA_WITH_THREADED_SOLVER
 #include "@DCA_THREADED_SOLVER_INCLUDE@"
-template<dca::DistType DIST>
+template <dca::DistType DIST>
 using ClusterSolver = @DCA_THREADED_SOLVER_TYPE@;
 #else
-template<dca::DistType DIST>
+template <dca::DistType DIST>
 using ClusterSolver = ClusterSolverBaseType<DIST>;
 #endif  // DCA_WITH_THREADED_SOLVER
 
 #include "dca/phys/dca_loop/dca_loop.hpp"
-template<dca::DistType DIST>
+template <dca::DistType DIST>
 using DcaLoopType = dca::phys::DcaLoop<ParametersType, DcaDataType<DIST>, ClusterSolver<DIST>, DIST>;
 
 #endif  // DCA_CONFIG_DCA_HPP

include/dca/config/mc_options.hpp.in

@@ -23,9 +23,9 @@ namespace config {
 // dca::config::
 
 struct McOptions {
-  using MCScalar = @MC_SCALAR@;
+  using MC_REAL = @MC_REAL@;
 
-  using TPAccumulationScalar = @TP_ACCUMULATION_SCALAR@;
+  using TPAccumulationPrecision = @TP_ACCUMULATION_PRECISION@;
 
   static constexpr bool memory_savings = @MEMORY_SAVINGS@;
 

include/dca/function/function.hpp

@@ -34,6 +34,7 @@
 #include "dca/util/pack_operations.hpp"
 #include "dca/util/integer_division.hpp"
 #include "dca/util/type_utils.hpp"
+#include "dca/util/type_help.hpp"
 #include "dca/util/to_string.hpp"
 
 namespace dca {
@@ -92,6 +93,9 @@ class function {
   // Same as above, but with name change from name argument.
   function(const function<scalartype, domain, DT>& other, const std::string& name);
 
+  template <typename Scalar2>
+  function(const function<Scalar2, domain, DT>& other);
+
   // Move constructor
   // 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
@@ -314,7 +318,6 @@ class function {
     return fnc_values_[dmn(static_cast<int>(t), static_cast<int>(subindices)...)];
   }
 
-
   void operator+=(const function<scalartype, domain, DT>& other);
   void operator-=(const function<scalartype, domain, DT>& other);
   void operator*=(const function<scalartype, domain, DT>& other);
@@ -363,7 +366,6 @@ class function {
   template <typename new_scalartype>
   void distribute(int sbdm_index_1, int sbdm_index_2, int* subind, const new_scalartype* fnc_vals);
 
-
   //
   // Methods for printing
   //
@@ -426,6 +428,13 @@ function<scalartype, domain, DT>::function(const std::string& name)
   }
   start_ = 0;
   end_ = dmn.get_size();
+  // If the function is more than 256 megs report it.
+  if (end_ > 268435456) {
+    std::cerr << "function " << name << " allocates " << sizeof(scalartype) * end_ / 1024 / 1024
+              << " MB" << '\n';
+    if (name_ == "no-name")
+      std::cerr << "large functions need names give yourself a chance.\n";
+  }
   // will zero real or complex values
   fnc_values_.resize(dmn.get_size(), {});
 }
@@ -457,6 +466,25 @@ function<scalartype, domain, DT>::function(const function<scalartype, domain, DT
   end_ = other.end_;
 }
 
+/** converting "copy" constructor
+ */
+template <typename scalartype, class domain, DistType DT>
+template <typename Scalar2>
+function<scalartype, domain, DT>::function(const function<Scalar2, domain, DT>& other)
+    : name_(other.get_name()),
+      function_type(__PRETTY_FUNCTION__),
+      dmn(),
+      Nb_sbdms(dmn.get_leaf_domain_sizes().size()),
+      fnc_values_(dmn.get_size()) {
+  if (size() != other.size()) {
+    // The other function has not been resetted after the domain was initialized.
+    throw std::logic_error("Copy construction from a not yet resetted function.");
+  }
+  start_ = other.get_start();
+  end_ = other.get_end();
+  std::copy(other.begin(), other.end(), begin());
+}
+
 /** move constructor */
 template <typename scalartype, class domain, DistType DT>
 function<scalartype, domain, DT>::function(function<scalartype, domain, DT>&& other)
@@ -640,7 +668,23 @@ template <typename Scalar, class domain, DistType DT>
 template <typename Scalar2>
 inline function<Scalar, domain, DT>& function<Scalar, domain, DT>::operator=(
     const function<Scalar2, domain, DT>& other) {
-  if (this != &other) {
+  if constexpr (std::is_same_v<decltype(*this), decltype(other)>) {
+    if (this != &other) {
+      if constexpr (dist == DistType::NONE) {
+        if (size() != other.size()) {
+          throw(std::logic_error("Function size does not match."));
+        }
+      }
+      else if constexpr (dist == DistType::LINEAR || dist == DistType::BLOCKED) {
+        Nb_sbdms = other.dmn.get_leaf_domain_sizes().size();
+        start_ = other.start_;
+        end_ = other.end_;
+        fnc_values_.resize(other.size(), {});
+      }
+      fnc_values_ = other.fnc_values_;
+    }
+  }
+  else {
     if constexpr (dist == DistType::NONE) {
       if (size() != other.size()) {
         throw(std::logic_error("Function size does not match."));
@@ -652,7 +696,16 @@ inline function<Scalar, domain, DT>& function<Scalar, domain, DT>::operator=(
       end_ = other.end_;
       fnc_values_.resize(other.size(), {});
     }
-    fnc_values_ = other.fnc_values_;
+    auto kConvert = [](auto& kvec) -> std::vector<Scalar> {
+      std::vector<Scalar> k_converted(kvec.size());
+      std::transform(kvec.begin(), kvec.end(), k_converted.begin(),
+                     [](auto& val) -> typename decltype(k_converted)::value_type {
+                       return static_cast<typename decltype(k_converted)::value_type>(val);
+                     });
+      return k_converted;
+    };
+
+    fnc_values_ = kConvert(other.getValues());
   }
   return *this;
 }

include/dca/io/adios2/adios2_reader.hpp

@@ -20,6 +20,7 @@
 
 #include "adios2.h"
 
+#include "dca/platform/dca_gpu.h"
 #include "dca/io/buffer.hpp"
 #include "dca/function/domains.hpp"
 #include "dca/function/function.hpp"