Add Vpd and Vpp to three-band Hubbard model

include/dca/phys/models/analytic_hamiltonians/threeband_hubbard.hpp

@@ -14,13 +14,15 @@
 #ifndef DCA_PHYS_MODELS_ANALYTIC_HAMILTONIANS_THREEBAND_LATTICE_HPP
 #define DCA_PHYS_MODELS_ANALYTIC_HAMILTONIANS_THREEBAND_LATTICE_HPP
 
+#include <algorithm>
 #include <cmath>
 #include <stdexcept>
 #include <utility>
 #include <vector>
 
 #include "dca/function/domains.hpp"
 #include "dca/function/function.hpp"
+#include "dca/phys/domains/cluster/cluster_operations.hpp"
 #include "dca/phys/domains/cluster/symmetries/point_groups/2d/2d_square.hpp"
 #include "dca/phys/domains/cluster/symmetries/point_groups/no_symmetry.hpp"
 #include "dca/phys/models/analytic_hamiltonians/cluster_shape_type.hpp"
@@ -31,7 +33,6 @@ namespace phys {
 namespace models {
 // dca::phys::models::
 
-// TODO: the symmetry of this system must be checked.
 template <typename SymmetryGroup>
 class ThreebandHubbard {
 public:
@@ -40,7 +41,6 @@ class ThreebandHubbard {
 
   using LDA_point_group = domains::no_symmetry<2>;
   using DCA_point_group = SymmetryGroup;
-  // typedef PointGroupType DCA_point_group;
 
   const static ClusterShapeType DCA_cluster_shape = BETT_CLUSTER;
   const static ClusterShapeType LDA_cluster_shape = PARALLELEPIPED;
@@ -59,11 +59,9 @@ class ThreebandHubbard {
 
   static std::vector<std::pair<std::pair<int, int>, std::pair<int, int>>> orbitalPermutations();
 
-  // Rotations of pi/2 are an anti-symmetry on the band off-diagonal.
   static int transformationSignOfR(int b1, int b2, int s);
   static int transformationSignOfK(int b1, int b2, int s);
 
-  // Initializes the interaction Hamiltonian in real space.
   template <typename BandDmn, typename SpinDmn, typename RDmn, typename parameters_type>
   static void initializeHInteraction(
       func::function<double, func::dmn_variadic<func::dmn_variadic<BandDmn, SpinDmn>,
@@ -72,9 +70,7 @@ class ThreebandHubbard {
 
   template <class domain>
   static void initializeHSymmetry(func::function<int, domain>& H_symmetry);
-
-  // Initializes the tight-binding (non-interacting) part of the momentum space Hamiltonian.
-  // Preconditions: The elements of KDmn are two-dimensional (access through index 0 and 1).
+
   template <typename ParametersType, typename ScalarType, typename BandDmn, typename SpinDmn, typename KDmn>
   static void initializeH0(
       const ParametersType& parameters,
@@ -110,7 +106,7 @@ int ThreebandHubbard<PointGroupType>::transformationSignOfR(int b1, int b2, int
       return 0;
   }
 
-  return s == 0;  // Only identity by default.
+  return s == 0;
 }
 
 template <typename PointGroupType>
@@ -120,19 +116,16 @@ int ThreebandHubbard<PointGroupType>::transformationSignOfK(int b1, int b2, int
 
   if ((b1 == b2) || (b1 != 0 && b2 != 0))
     return 1;
-  else  // if ((b1 != b2 && b1 == 0)||(b1 != b2 && b2 == 0))
+  else
     return (s == 0 || s == 2 || s == 5 || s == 7) ? 1 : -1;
-  // if (s == 0 || s == 2 || s == 5 || s == 7)
-  //  return 1;
-  // else // if (s == 1 || s == 3 || s == 4 || s == 6)
-  //  return -1;
 }
 
 template <typename PointGroupType>
 const double* ThreebandHubbard<PointGroupType>::initializeRDCABasis() {
   static std::array<double, 4> basis{1, 0, 0, 1};
   return basis.data();
 }
+
 template <typename PointGroupType>
 const double* ThreebandHubbard<PointGroupType>::initializeKDCABasis() {
   static std::array<double, 4> basis{2 * M_PI, 0, 0, 2 * M_PI};
@@ -162,8 +155,8 @@ std::vector<std::vector<double>> ThreebandHubbard<PointGroupType>::aVectors() {
 }
 
 template <typename PointGroupType>
-std::vector<std::pair<std::pair<int, int>, std::pair<int, int>>> ThreebandHubbard<
-    PointGroupType>::orbitalPermutations() {
+std::vector<std::pair<std::pair<int, int>, std::pair<int, int>>>
+ThreebandHubbard<PointGroupType>::orbitalPermutations() {
   return {};
 }
 
@@ -180,11 +173,65 @@ void ThreebandHubbard<PointGroupType>::initializeHInteraction(
 
   const int origin = RDmn::parameter_type::origin_index();
 
-  const double U_dd = parameters.get_U_dd();  // interaction in d band
-  const double U_pp = parameters.get_U_pp();  // interaction in p bands
+  const double U_dd = parameters.get_U_dd();
+  const double U_pp = parameters.get_U_pp();
+  const double V_pd = parameters.get_V_pd();
+  const double V_pp = parameters.get_V_pp();
 
   H_interaction = 0.;
 
+  const auto& basis = RDmn::parameter_type::get_basis_vectors();
+  const auto& super_basis = RDmn::parameter_type::get_super_basis_vectors();
+  const auto& elements = RDmn::parameter_type::get_elements();
+
+  if (basis.size() != DIMENSION)
+    throw std::logic_error("Unexpected lattice basis dimension for three-band Hubbard model.");
+
+  auto index_of = [&](std::vector<double> displacement) {
+    displacement = domains::cluster_operations::translate_inside_cluster(displacement, super_basis);
+    return domains::cluster_operations::index(displacement, elements, domains::BRILLOUIN_ZONE);
+  };
+
+  auto push_unique = [](std::vector<int>& indices, int index) {
+    if (std::find(indices.begin(), indices.end(), index) == indices.end())
+      indices.push_back(index);
+  };
+
+  const int plus_x = index_of(basis[0]);
+  const int plus_y = index_of(basis[1]);
+  const int minus_x = RDmn::parameter_type::subtract(plus_x, origin);
+  const int minus_y = RDmn::parameter_type::subtract(plus_y, origin);
+
+  std::vector<double> plus_x_minus_y(DIMENSION, 0.);
+  std::vector<double> minus_x_plus_y(DIMENSION, 0.);
+  for (int d = 0; d < DIMENSION; ++d) {
+    plus_x_minus_y[d] = basis[0][d] - basis[1][d];
+    minus_x_plus_y[d] = basis[1][d] - basis[0][d];
+  }
+
+  const int delta_plus_x_minus_y = index_of(plus_x_minus_y);
+  const int delta_minus_x_plus_y = index_of(minus_x_plus_y);
+
+  std::vector<int> d_to_px{origin};
+  std::vector<int> px_to_d{origin};
+  std::vector<int> d_to_py{origin};
+  std::vector<int> py_to_d{origin};
+  std::vector<int> px_to_py{origin};
+  std::vector<int> py_to_px{origin};
+
+  push_unique(d_to_px, plus_x);
+  push_unique(px_to_d, minus_x);
+  push_unique(d_to_py, plus_y);
+  push_unique(py_to_d, minus_y);
+
+  push_unique(px_to_py, minus_x);
+  push_unique(px_to_py, plus_y);
+  push_unique(px_to_py, delta_minus_x_plus_y);
+
+  push_unique(py_to_px, plus_x);
+  push_unique(py_to_px, minus_y);
+  push_unique(py_to_px, delta_plus_x_minus_y);
+
   for (int i = 0; i < BANDS; i++) {
     for (int s1 = 0; s1 < 2; s1++) {
       for (int j = 0; j < BANDS; j++) {
@@ -198,37 +245,47 @@ void ThreebandHubbard<PointGroupType>::initializeHInteraction(
       }
     }
   }
+
+  for (int s1 = 0; s1 < 2; s1++) {
+    for (int s2 = 0; s2 < 2; s2++) {
+      for (const int delta_r : d_to_px)
+        H_interaction(0, s1, 1, s2, delta_r) = V_pd;
+      for (const int delta_r : d_to_py)
+        H_interaction(0, s1, 2, s2, delta_r) = V_pd;
+      for (const int delta_r : px_to_d)
+        H_interaction(1, s1, 0, s2, delta_r) = V_pd;
+      for (const int delta_r : py_to_d)
+        H_interaction(2, s1, 0, s2, delta_r) = V_pd;
+    }
+  }
+
+  for (int s1 = 0; s1 < 2; s1++) {
+    for (int s2 = 0; s2 < 2; s2++) {
+      for (const int delta_r : px_to_py)
+        H_interaction(1, s1, 2, s2, delta_r) = V_pp;
+      for (const int delta_r : py_to_px)
+        H_interaction(2, s1, 1, s2, delta_r) = V_pp;
+    }
+  }
 }
 
 template <typename PointGroupType>
 template <class domain>
-void ThreebandHubbard<PointGroupType>::initializeHSymmetry(func::function<int, domain>& H_symmetries) {
+void ThreebandHubbard<PointGroupType>::initializeHSymmetry(
+    func::function<int, domain>& H_symmetries) {
   H_symmetries = -1;
-
-  //  H_symmetry(i, s1, j, s2)
-  //  H_symmetries(0, 0, 0, 0) = 0; // at b0, G of spin 0 or 1 has the same values.
-  //  H_symmetries(0, 1, 0, 1) = 0;
-
-  //  H_symmetries(1, 0, 1, 0) = 1;
-  //  H_symmetries(1, 1, 1, 1) = 1; // at i, G of spin 0 or 1 has the same values.
 }
 
-
 template <typename PointGroupType>
 template <typename ParametersType, typename ScalarType, typename BandDmn, typename SpinDmn, typename KDmn>
 void ThreebandHubbard<PointGroupType>::initializeH0(
     const ParametersType& parameters,
     func::function<ScalarType, func::dmn_variadic<func::dmn_variadic<BandDmn, SpinDmn>,
     func::dmn_variadic<BandDmn, SpinDmn>, KDmn>>& H_0) {
   typename KDmn::element_type default_q(ParametersType::lattice_dimension);
-  // would rather get it like this -> KDmn::parameter_type::DIMENSION);
-  // but KDmn are inconsistent about whether they known their kspace dimension!
   initializeH0WithQ(parameters, H_0, default_q);
 }
 
-/** generate the actual H0 with q-shift
- *  assumption that the origin of KDmn is the default initialization of KDmn::element_type
- */
 template <typename PointGroupType>
 template <typename ParametersType, typename ScalarType, typename BandDmn, typename SpinDmn, typename KDmn>
 void ThreebandHubbard<PointGroupType>::initializeH0WithQ(
@@ -242,29 +299,34 @@ void ThreebandHubbard<PointGroupType>::initializeH0WithQ(
 
   std::vector<typename KDmn::element_type> k_vecs(KDmn::get_elements());
 
-  for( auto & k_elem : k_vecs)
-    std::transform(k_elem.cbegin(), k_elem.cend(), q.begin(), k_elem.begin(), [](auto& k, auto q_elem) { return k + q_elem; });
-
+  for (auto& k_elem : k_vecs)
+    std::transform(k_elem.cbegin(), k_elem.cend(), q.begin(), k_elem.begin(),
+                   [](auto& k, auto q_elem) { return k + q_elem; });
+
   const auto t_pd = parameters.get_t_pd();
   const auto t_pp = parameters.get_t_pp();
   const auto ep_d = parameters.get_ep_d();
   const auto ep_p = parameters.get_ep_p();
+  const auto U_dd = parameters.get_U_dd();
+  const auto U_pp = parameters.get_U_pp();
+  const auto V_pd = parameters.get_V_pd();
+  const auto V_pp = parameters.get_V_pp();
 
   H_0 = ScalarType(0);
 
   const ScalarType I(0, 1.);
 
   for (int k_ind = 0; k_ind < KDmn::dmn_size(); ++k_ind) {
     const auto& k = k_vecs[k_ind];
-    
-    const auto valdpx = 2. * I * t_pd * std::sin(k[0] / 2.);
-    const auto valdpy = -2. * I * t_pd * std::sin(k[1] / 2.);
+
+    const auto valdpx = -2. * I * t_pd * std::sin(k[0] / 2.);
+    const auto valdpy = 2. * I * t_pd * std::sin(k[1] / 2.);
     const auto valpxpy = 4. * t_pp * std::sin(k[0] / 2.) * std::sin(k[1] / 2.);
 
     for (int s = 0; s < 2; s++) {
-      H_0(0, s, 0, s, k_ind) = ep_d;
-      H_0(1, s, 1, s, k_ind) = ep_p;
-      H_0(2, s, 2, s, k_ind) = ep_p;
+      H_0(0, s, 0, s, k_ind) = ep_d + U_dd / 2.0 + 4.0 * V_pd;
+      H_0(1, s, 1, s, k_ind) = ep_p + U_pp / 2.0 + 4.0 * V_pp + 2.0 * V_pd;
+      H_0(2, s, 2, s, k_ind) = ep_p + U_pp / 2.0 + 4.0 * V_pp + 2.0 * V_pd;
 
       H_0(0, s, 1, s, k_ind) = valdpx;
       H_0(1, s, 0, s, k_ind) = -valdpx;

include/dca/phys/parameters/model_parameters_threeband_hubbard.inc

@@ -67,13 +67,29 @@ public:
     U_pp_ = U_pp;
   }
 
+  double get_V_pd() const {
+    return V_pd_;
+  }
+  void set_V_pd(const double V_pd) {
+    V_pd_ = V_pd;
+  }
+
+  double get_V_pp() const {
+    return V_pp_;
+  }
+  void set_V_pp(const double V_pp) {
+    V_pp_ = V_pp;
+  }
+
 private:
   double t_pd_ = 0;
   double t_pp_ = 0;
   double ep_p_ = 0;
   double ep_d_ = 0;
   double U_dd_ = 0;
   double U_pp_ = 0;
+  double V_pd_ = 0;
+  double V_pp_ = 0;
 };
 
 template <typename PointGroup>
@@ -88,6 +104,8 @@ int ModelParameters<models::TightBindingModel<models::ThreebandHubbard<PointGrou
   buffer_size += concurrency.get_buffer_size(ep_d_);
   buffer_size += concurrency.get_buffer_size(U_dd_);
   buffer_size += concurrency.get_buffer_size(U_pp_);
+  buffer_size += concurrency.get_buffer_size(V_pd_);
+  buffer_size += concurrency.get_buffer_size(V_pp_);
 
   return buffer_size;
 }
@@ -102,6 +120,8 @@ void ModelParameters<models::TightBindingModel<models::ThreebandHubbard<PointGro
   concurrency.pack(buffer, buffer_size, position, ep_d_);
   concurrency.pack(buffer, buffer_size, position, U_dd_);
   concurrency.pack(buffer, buffer_size, position, U_pp_);
+  concurrency.pack(buffer, buffer_size, position, V_pd_);
+  concurrency.pack(buffer, buffer_size, position, V_pp_);
 }
 
 template <typename PointGroup>
@@ -114,6 +134,8 @@ void ModelParameters<models::TightBindingModel<models::ThreebandHubbard<PointGro
   concurrency.unpack(buffer, buffer_size, position, ep_d_);
   concurrency.unpack(buffer, buffer_size, position, U_dd_);
   concurrency.unpack(buffer, buffer_size, position, U_pp_);
+  concurrency.unpack(buffer, buffer_size, position, V_pd_);
+  concurrency.unpack(buffer, buffer_size, position, V_pp_);
 }
 
 template <typename PointGroup>
@@ -153,6 +175,16 @@ void ModelParameters<models::TightBindingModel<models::ThreebandHubbard<PointGro
     }
     catch (const std::exception& r_e) {
     }
+    try {
+      reader_or_writer.execute("V_pd", V_pd_);
+    }
+    catch (const std::exception& r_e) {
+    }
+    try {
+      reader_or_writer.execute("V_pp", V_pp_);
+    }
+    catch (const std::exception& r_e) {
+    }
 
     reader_or_writer.close_group();
   }

test/unit/phys/models/analytic_hamiltonians/CMakeLists.txt

@@ -19,3 +19,8 @@ dca_add_gtest(triangular_lattice_test
   GTEST_MAIN
   INCLUDE_DIRS ${FFTW_INCLUDE_DIR}
   LIBS function function_transform ${LAPACK_LIBRARIES} ${DCA_GPU_LIBS}  ${UTIL_LIBS})
+
+dca_add_gtest(threeband_hubbard_test
+  GTEST_MAIN
+  INCLUDE_DIRS ${FFTW_INCLUDE_DIR}
+  LIBS function function_transform json ${LAPACK_LIBRARIES} ${DCA_GPU_LIBS} ${UTIL_LIBS})

test/unit/phys/models/analytic_hamiltonians/threeband_hubbard_input.json

@@ -0,0 +1,12 @@
+{
+  "threebands-Hubbard-model": {
+    "t_pd": 1.5,
+    "t_pp": 0.25,
+    "ep_d": -1.0,
+    "ep_p": 2.0,
+    "U_dd": 8.0,
+    "U_pp": 4.0,
+    "V_pd": 0.75,
+    "V_pp": 0.5
+  }
+}