merge

examples/dynamics_finite_t/shastry_sutherland_fulled/plot_spectra.py

@@ -63,9 +63,6 @@
     poles = poles[np.abs(weights) > min_weight]
     weights = np.real(weights[np.abs(weights) > min_weight])
 
-    print(weights)
-    print(np.sort(poles))
-
     # broaden and plot
     diffs = np.subtract.outer(omegas, poles)
     gaussians = np.exp(-(diffs / (2*eta))**2) / (eta * np.sqrt(2*np.pi))

hydra/CMakeLists.txt

@@ -15,7 +15,9 @@ add_library(hydra STATIC
 
   algebra/algebra.cpp
   algebra/matrix.cpp
-
+
+  io/args.cpp
+  io/args_handler.cpp
   io/file_toml.cpp
   io/toml/file_toml_handler.cpp
   io/toml/toml_conversion.cpp
@@ -71,7 +73,8 @@ add_library(hydra STATIC
   blocks/tj/tj_apply.cpp
   blocks/tj/terms/compile.cpp
 
-  symmetries/operations/symmetry_operations.cpp
+  symmetries/qn.cpp
+  symmetries/operations/symmetry_operations.cpp	
   symmetries/permutation.cpp
   symmetries/permutation_group.cpp
   symmetries/generated_group.cpp

hydra/all.h

@@ -89,6 +89,8 @@
 #include "blocks/tj/tj_matrix.h"
 #include "blocks/tj/tj_utils.h"
 
+#include "symmetries/qn.h"
+#include "symmetries/continuous_group.h"
 #include "symmetries/generated_group.h"
 #include "symmetries/group_action/group_action.h"
 #include "symmetries/group_action/group_action_lookup.h"
@@ -128,6 +130,7 @@
 
 #include "io/file_h5.h"
 #include "io/file_toml.h"
+#include "io/args.h"
 
 #ifdef HYDRA_ENABLE_MPI
 #include "mpi/allreduce.h"

hydra/blocks/spinhalf/spinhalf_matrix.cpp

@@ -7,6 +7,7 @@
 #include <hydra/blocks/utils/block_utils.h>
 #include <hydra/operators/compiler.h>
 #include <hydra/utils/logger.h>
+#include <hydra/utils/timing.h>
 
 namespace hydra {
 
@@ -22,13 +23,13 @@ arma::Mat<coeff_t> matrix_gen(BondList const &bonds, Spinhalf const &block_in,
             "intrisically complex BondList");
   }
 
-  if (block_in.n_up() != undefined_qn) {
-    int n_up_out = spinhalf::nup(bonds_c) + block_in.n_up();
-    if (n_up_out != block_out.n_up()) {
-      Log.err("Incompatible n_up in matrix_gen: {} != {}", n_up_out,
-              block_out.n_up());
-    }
-  }
+  // if (block_in.n_up() != undefined_qn) {
+  //   int n_up_out = spinhalf::nup(bonds_c) + block_in.n_up();
+  //   if (n_up_out != block_out.n_up()) {
+  //     Log.err("Incompatible n_up in matrix_gen: {} != {}", n_up_out,
+  //             block_out.n_up());
+  //   }
+  // }
 
   idx_t dim_in = block_in.size();
   idx_t dim_out = block_out.size();
@@ -42,7 +43,7 @@ arma::Mat<coeff_t> matrix_gen(BondList const &bonds, Spinhalf const &block_in,
   auto const &indexing_out = block_out.indexing();
   spinhalf::apply_terms_dispatch<coeff_t>(bonds_c, indexing_in, indexing_out,
                                           fill);
-  return mat;
+ return mat;
 }
 
 template arma::Mat<double> matrix_gen<double>(BondList const &bonds,

hydra/blocks/spinhalf/terms/apply_terms.h

@@ -10,6 +10,7 @@
 #include <hydra/common.h>
 
 #include <hydra/utils/print_macro.h>
+#include <hydra/utils/timing.h>
 
 namespace hydra::spinhalf {
 

hydra/blocks/spinhalf/terms/compile.cpp

@@ -3,6 +3,7 @@
 #include <hydra/operators/non_branching_bonds.h>
 #include <hydra/utils/logger.h>
 #include <hydra/utils/print_macro.h>
+#include <hydra/utils/timing.h>
 
 namespace hydra::spinhalf {
 

hydra/indexing/fermi_table.cpp

@@ -61,12 +61,10 @@ std::vector<bool> init_fermi_table_omp(States const &states,
       auto states_thread = ThreadStates(states);
       // fermi_bool_table_local[myid].resize(states_thread.size());
       auto fermi_work = symmetries::fermi_work(n_sites);
-      idx_t idx = 0;
       for (auto state : states_thread) {
         bool fermi_bool =
             symmetries::fermi_bool_of_permutation(state, perm, fermi_work);
         fermi_bool_table_local[myid].push_back(fermi_bool);
-        ++idx;
       }
 #pragma omp barrier
 

hydra/io/args.cpp

@@ -0,0 +1,19 @@
+#include "args.h"
+
+namespace hydra {
+
+Args::Args(
+    std::initializer_list<std::pair<std::string, io::args_t>> const &args) {
+  for (auto const &[key, val] : args) {
+    args_[key] = val;
+  }
+}
+bool Args::defined(std::string key) const { return args_.count(key); }
+io::ArgsHandler Args::operator[](std::string key) {
+  return io::ArgsHandler(key, args_);
+}
+
+bool Args::operator==(Args const &other) const { return args_ == other.args_; }
+bool Args::operator!=(Args const &other) const { return !operator==(other); }
+
+} // namespace hydra

hydra/io/args.h

@@ -0,0 +1,27 @@
+#pragma once
+
+#include <string>
+#include <map>
+#include <initializer_list>
+#include <utility>
+
+#include <hydra/io/args_handler.h>
+
+namespace hydra {
+
+class Args {
+public:
+  Args() = default;
+  Args(std::initializer_list<std::pair<std::string, io::args_t>> const& args);
+
+  bool defined(std::string key) const;
+  io::ArgsHandler operator[](std::string key);
+
+  bool operator==(Args const &other) const;
+  bool operator!=(Args const &other) const;
+
+private:
+  std::map<std::string, io::args_t> args_;
+};
+
+} // namespace hydra

hydra/io/args_handler.cpp

@@ -0,0 +1,70 @@
+#include "args_handler.h"
+
+#include <hydra/utils/logger.h>
+
+namespace hydra::io {
+
+ArgsHandler::ArgsHandler(std::string key,
+                         std::map<std::string, io::args_t> &args)
+    : key_(key), args_(args) {}
+
+template <class data_t> data_t ArgsHandler::as() const {
+  if (args_.count(key_)) {
+    return std::get<data_t>(args_.at(key_));
+  } else {
+    Log.err("Error using Args: key \"{}\" not found", key_);
+    return data_t();
+  }
+}
+template <class data_t> data_t ArgsHandler::as(data_t const &defaultt) const {
+  if (args_.count(key_)) {
+    return std::get<data_t>(args_.at(key_));
+  } else {
+    return defaultt;
+  }
+}
+
+template <class data_t> void ArgsHandler::operator=(data_t const &data) {
+  args_[key_] = data;
+}
+
+template bool ArgsHandler::as<bool>() const;
+template std::string ArgsHandler::as<std::string>() const;
+template int8_t ArgsHandler::as<int8_t>() const;
+template int16_t ArgsHandler::as<int16_t>() const;
+template int32_t ArgsHandler::as<int32_t>() const;
+template int64_t ArgsHandler::as<int64_t>() const;
+template uint8_t ArgsHandler::as<uint8_t>() const;
+template uint16_t ArgsHandler::as<uint16_t>() const;
+template uint32_t ArgsHandler::as<uint32_t>() const;
+template uint64_t ArgsHandler::as<uint64_t>() const;
+template double ArgsHandler::as<double>() const;
+template complex ArgsHandler::as<complex>() const;
+
+template bool ArgsHandler::as(bool const &) const;
+template std::string ArgsHandler::as(std::string const &) const;
+template int8_t ArgsHandler::as(int8_t const &) const;
+template int16_t ArgsHandler::as(int16_t const &) const;
+template int32_t ArgsHandler::as(int32_t const &) const;
+template int64_t ArgsHandler::as(int64_t const &) const;
+template uint8_t ArgsHandler::as(uint8_t const &) const;
+template uint16_t ArgsHandler::as(uint16_t const &) const;
+template uint32_t ArgsHandler::as(uint32_t const &) const;
+template uint64_t ArgsHandler::as(uint64_t const &) const;
+template double ArgsHandler::as(double const &) const;
+template complex ArgsHandler::as(complex const &) const;
+
+template void ArgsHandler::operator=(bool const &);
+template void ArgsHandler::operator=(std::string const &);
+template void ArgsHandler::operator=(int8_t const &);
+template void ArgsHandler::operator=(int16_t const &);
+template void ArgsHandler::operator=(int32_t const &);
+template void ArgsHandler::operator=(int64_t const &);
+template void ArgsHandler::operator=(uint8_t const &);
+template void ArgsHandler::operator=(uint16_t const &);
+template void ArgsHandler::operator=(uint32_t const &);
+template void ArgsHandler::operator=(uint64_t const &);
+template void ArgsHandler::operator=(double const &);
+template void ArgsHandler::operator=(complex const &);
+
+} // namespace hydra::io

hydra/io/args_handler.h

@@ -0,0 +1,32 @@
+#pragma once
+
+#include <map>
+#include <string>
+#include <variant>
+#include <vector>
+
+#include <extern/armadillo/armadillo>
+#include <hydra/common.h>
+
+namespace hydra::io {
+
+using args_t =
+    std::variant<int8_t, int16_t, int32_t, int64_t, uint8_t, uint16_t, uint32_t,
+                 uint64_t, bool, double, complex, std::string>;
+
+class ArgsHandler {
+public:
+  ArgsHandler(std::string key, std::map<std::string, io::args_t> &args);
+  ArgsHandler(ArgsHandler const &) = delete;
+  ArgsHandler &operator=(ArgsHandler const &) = delete;
+
+  template <class data_t> data_t as() const;
+  template <class data_t> data_t as(data_t const& defaultt) const;
+  template <class data_t> void operator=(data_t const &data);
+
+private:
+  std::string key_;
+  std::map<std::string, io::args_t> &args_;
+};
+
+} // namespace hydra::io

hydra/operators/non_branching_bonds.cpp

@@ -149,7 +149,7 @@ NonBranchingBond<bit_t, coeff_t>::NonBranchingBond(Bond const &bond,
   coeff_ = std::vector<coeff_t>(dim_, 0.);
 
   for (bit_t in = 0; in < dim_; ++in) {
-    int non_zero_in_row = 0;
+    // int non_zero_in_row = 0;
 
     for (bit_t out = 0; out < dim_; ++out) {
       if (std::abs(matrix_(out, in)) > precision) {
@@ -164,16 +164,16 @@ NonBranchingBond<bit_t, coeff_t>::NonBranchingBond(Bond const &bond,
         } else {
           coeff_[in] = cpl * matrix_(out, in);
         }
-        ++non_zero_in_row;
+        // ++non_zero_in_row;
       }
     }
 
-    // security check
-    if (non_zero_term_[in]) {
-      assert(non_zero_in_row == 1);
-    } else {
-      assert(non_zero_in_row == 0);
-    }
+    // // security check
+    // if (non_zero_term_[in]) {
+    //   assert(non_zero_in_row == 1);
+    // } else {
+    //   assert(non_zero_in_row == 0);
+    // }
   }
 
   // int n_sites = 1;

hydra/states/product_state.cpp

@@ -13,7 +13,7 @@ void fill(ProductState const &pstate, State<coeff_t> &state) {
   std::visit(variant::overloaded{[&pstate, &vector](auto &&block) {
                fill(pstate, block, vector);
              }},
-    state.block().variant());
+             state.block().variant());
 }
 template void fill(ProductState const &pstate, State<double> &state);
 template void fill(ProductState const &pstate, State<complex> &state);
@@ -142,7 +142,7 @@ std::pair<bit_t, bit_t> get_tj_spins(ProductState const &pstate, int nup = -1,
     }
   }
 
-  if ((nup != pnup) && (nup != -1) && (ndn != -1)) {
+  if (((nup != pnup) || (ndn != pndn)) && (nup != -1) && (ndn != -1)) {
     Log.err("Error creating product state: number of up/dn spins incompatible "
             "with block");
   }
@@ -218,7 +218,7 @@ std::pair<bit_t, bit_t> get_electron_spins(ProductState const &pstate,
     }
   }
 
-  if ((nup != pnup) && (nup != -1) && (ndn != -1)) {
+  if (((nup != pnup) || (ndn != pndn)) && (nup != -1) && (ndn != -1)) {
     Log.err("Error creating product state: number of up/dn spins incompatible "
             "with block");
   }

hydra/symmetries/continuous_group.h

@@ -0,0 +1,14 @@
+#pragma once
+
+namespace hydra {
+
+class U1 {
+public:
+  inline bool operator==(U1 const &rhs) const {
+    (void)rhs;
+    return true;
+  }
+  inline bool operator!=(U1 const &rhs) const { return !operator==(rhs); }
+};
+
+} // namespace hydra