Model Hamiltonians

include/macis/model/hubbard.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+#include <macis/types.hpp>
+
+namespace macis {
+
+/** 
+ *  @brief Generate Hamiltonian Data for 1D Hubbard
+ */
+void hubbard_1d(size_t nsites, double t, double U, std::vector<double>& T, 
+                std::vector<double>& V);
+
+}

src/macis/CMakeLists.txt

@@ -15,6 +15,7 @@ add_library( macis
   orbital_rotation_utilities.cxx
   orbital_hessian.cxx
   orbital_steps.cxx
+  model/hubbard.cxx
 )
 
 target_include_directories( macis PUBLIC 

src/macis/model/hubbard.cxx

@@ -0,0 +1,29 @@
+#include <macis/model/hubbard.hpp>
+
+namespace macis {
+
+void hubbard_1d(size_t nsites, double t, double U, std::vector<double>& T, 
+                std::vector<double>& V) {
+
+
+  T.resize(nsites * nsites);
+  V.resize(nsites * nsites * nsites * nsites);
+
+
+  for(size_t p = 0; p < nsites; ++p) {
+    // Half-filling Chemical Potential
+    T[p * (nsites+1)] = -U/2; 
+
+    // On-Site Interaction
+    V[p * (nsites*nsites*nsites + nsites*nsites + nsites + 1)] = U; 
+
+    // Hopping
+    if(p < nsites-1) {
+      T[p + (p+1)*nsites] = -t;
+      T[(p+1) + p*nsites] = -t;
+    }
+  }
+
+}
+
+}

tests/CMakeLists.txt

@@ -20,6 +20,7 @@ add_executable( macis_test
   mcscf.cxx
   asci.cxx
   dist_quickselect.cxx
+  hubbard.cxx
 )
 target_link_libraries( macis_test PUBLIC macis Catch2::Catch2 )
 target_include_directories( macis_test PUBLIC ${PROJECT_BINARY_DIR}/tests )

tests/hubbard.cxx

@@ -0,0 +1,48 @@
+
+
+#include <iostream>
+#include <spdlog/sinks/null_sink.h>
+#include <spdlog/spdlog.h>
+
+#include <macis/model/hubbard.hpp>
+
+#include "ut_common.hpp"
+
+TEST_CASE("Hubbard") {
+  ROOT_ONLY(MPI_COMM_WORLD);
+
+  const size_t nsites = 4;
+  const size_t nsites2 = nsites*nsites;
+  const size_t nsites3 = nsites2*nsites;
+  const double t = 1.0, U = 4.0;
+  std::vector<double> T,V;
+
+  SECTION("1D") {
+    macis::hubbard_1d( nsites, t, U, T, V );
+
+    // Check two-body term
+    for(int p = 0; p < nsites; ++p)
+    for(int q = 0; q < nsites; ++q)
+    for(int r = 0; r < nsites; ++r)
+    for(int s = 0; s < nsites; ++s) {
+      const auto mat_el = V[p + nsites*q + nsites2*r + nsites3*s];
+      if(p==q and p==r and p==s)
+        REQUIRE(mat_el == U);
+      else
+        REQUIRE(mat_el == 0.0);
+    }
+
+    // Check 1-body term
+    for(int p = 0; p < nsites; ++p)
+    for(int q = 0; q < nsites; ++q) {
+      const auto mat_el = T[p + q*nsites];
+      if(p == q) 
+        REQUIRE(mat_el == Approx(-U/2));
+      else if( std::abs(p-q) == 1 )
+        REQUIRE(mat_el == -t);
+      else
+        REQUIRE(mat_el == 0.0);
+    }
+  }
+
+}