Enable D0_iw and Jperp_iw terms

c++/triqs_ctint/solver_core.cpp

@@ -169,8 +169,8 @@ namespace triqs_ctint {
     g_iw_t G0_shift_iw_inv = G0_iw_inv;
 
     // Assert compatibility between alpha tensor and h_int
-    long n_terms = std::distance(p.h_int.begin(), p.h_int.end());
-    if (p.alpha.shape() != std::array<long, 4>{n_terms, 2, 2, p.n_s})
+    long const n_terms = std::distance(p.h_int.begin(), p.h_int.end());
+    if (p.alpha.shape() != std::array<long, 4>{n_terms + (D0_iw ? 1 : 0) + (Jperp_iw ? 1 : 0), 2, 2, p.n_s})
       TRIQS_RUNTIME_ERROR << "Error: Alpha and h_int are incompatible: Different number of density-density terms \n";
 
     // Loop over static density-density interaction terms
@@ -198,21 +198,19 @@ namespace triqs_ctint {
 
     if (D0_iw) {
 
-      // External loop over blocks
-      for (int sig : range(p.n_blocks())) {
-
-        // Get Matrix Rank for block
-        int Rank = G0_iw[sig].target_shape()[0];
-        for (int i : range(Rank)) {
-
-          // Calculate and subtract term according to notes
-          TRIQS_RUNTIME_ERROR << "FIXME";
-          //dcomplex term = 0.0;
-          //for (int sigp : range(p.n_blocks()))
-          //for (int j : range(Rank))
-          //for (int s : range(p.n_s)) { term += ((*D0_iw)(sig, sigp)[0](i, j) + (*D0_iw)(sigp, sig)[0](j, i)) * p.alpha_l[sigp](j, s); }
-          //auto g = slice_target_to_scalar(G0_shift_iw_inv[sig], i, i);
-          //g(iw_) << g(iw_) - term / p.n_s;
+      for (int bl1 : range(p.n_blocks())) {
+        int bl1_size = G0_iw[bl1].target_shape()[0];
+        for (int i : range(bl1_size)) {
+          dcomplex term = 0.0;
+          for (int bl2 : range(p.n_blocks())) {
+            int bl2_size = G0_iw[bl2].target_shape()[0];
+            for (int j : range(bl2_size)) {
+              for (int s : range(p.n_s)) {
+                term += ((*D0_iw)(bl1, bl2)[0](i, j) + (*D0_iw)(bl2, bl1)[0](j, i)) * p.alpha(n_terms, bl2, bl2, s);
+              }
+            }
+          }
+          G0_shift_iw_inv[bl1].data()(range(), i, i) -= term / p.n_s;
         }
       }
     }

c++/triqs_ctint/vertex_factories.cpp

@@ -9,6 +9,10 @@ namespace triqs_ctint {
 
     std::vector<vertex_factory_t> vertex_factories;
 
+    int const n_terms = std::distance(params.h_int.begin(), params.h_int.end());
+    bool const use_D = static_cast<bool>(D0_iw);
+    bool const use_Jperp = static_cast<bool>(Jperp_iw);
+
     // ------------ Create Vertex Factory for Static Interactions --------------
     {
       std::vector<vertex_idx_t> indices;
@@ -44,7 +48,7 @@ namespace triqs_ctint {
     } // clean temporaries
 
     // ------------ Create Vertex Factory for Dynamic Density-Density Interactions --------------
-    if (D0_iw) {
+    if (use_D) {
 
       std::vector<vertex_idx_t> indices;
 #ifdef INTERACTION_IS_COMPLEX
@@ -80,22 +84,22 @@ namespace triqs_ctint {
         }
 
       if (indices.size() > 0) {
-        auto l = [beta = params.beta, n_s = params.n_s, indices = std::move(indices), D0_tau_lst = std::move(D0_tau_lst), &rng] {
+        auto l = [beta = params.beta, n_s = params.n_s, indices = std::move(indices), D0_tau_lst = std::move(D0_tau_lst), &rng, n_terms] {
           int n             = rng(indices.size());
           tau_t t           = tau_t::get_random(rng);
           tau_t tp          = tau_t::get_random(rng);
           auto [sig, dtau]  = cyclic_difference(t, tp);
           int s             = rng(n_s);
           double prop_proba = 1.0 / (beta * beta * indices.size() * n_s);
-          return vertex_t{indices[n], t, t, tp, tp, -D0_tau_lst[n](dtau) / n_s, prop_proba, true, s};
+          return vertex_t{indices[n], t, t, tp, tp, -D0_tau_lst[n](dtau) / n_s, prop_proba, n_terms, s};
         };
 
         vertex_factories.emplace_back(l);
       }
     }
 
     // ------------ Create Vertex Factory for Dynamic Spin-Spin Interactions --------------
-    if (Jperp_iw) {
+    if (use_Jperp) {
 
       std::vector<vertex_idx_t> indices;
 #ifdef INTERACTION_IS_COMPLEX
@@ -129,13 +133,13 @@ namespace triqs_ctint {
         }
 
       if (indices.size() > 0) {
-        auto l = [beta = params.beta, indices = std::move(indices), Jperp_tau_lst = std::move(Jperp_tau_lst), &rng] {
+        auto l = [beta = params.beta, indices = std::move(indices), Jperp_tau_lst = std::move(Jperp_tau_lst), &rng, n_terms, use_D] {
           int n             = rng(indices.size());
           tau_t t           = tau_t::get_random(rng);
           tau_t tp          = tau_t::get_random(rng);
           auto [sig, dtau]  = cyclic_difference(t, tp);
           double prop_proba = 1.0 / (beta * beta * indices.size());
-          return vertex_t{indices[n], t, t, tp, tp, Jperp_tau_lst[n](dtau) / 2.0, prop_proba}; // We add two identical terms above -> Divide by 2
+          return vertex_t{indices[n], t, t, tp, tp, Jperp_tau_lst[n](dtau) / 2.0, prop_proba, n_terms + use_D}; // We add two identical terms above -> Divide by 2
         };
 
         vertex_factories.emplace_back(l);

python/triqs_ctint/solver.py

@@ -172,6 +172,8 @@ def jacobi(self, x, solve_params):
     def find_alpha_from_self_consistent_HF(self, solve_params):
         # --------- Determine the alpha tensor from SC Hartree Fock ----------
         mpi_print("Determine alpha-tensor")
+        assert not (self.constr_params['use_D'] or self.constr_params['use_Jperp']), \
+            "Automatic determination of alpha tensor does not work with D0_iw or Jperp_iw"
 
         gf_struct = self.constr_params['gf_struct']
         h_int = solve_params['h_int']
@@ -273,12 +275,14 @@ def sign(a):
 
     def trivial_alpha(self, solve_params):
         gf_struct = self.constr_params['gf_struct']
+        use_D = self.constr_params['use_D']
+        use_Jperp = self.constr_params['use_Jperp']
         h_int = solve_params['h_int']
         n_terms = len(list(h_int))
         delta = solve_params.pop('delta', [0.5 + 1e-2, 1e-2])
 
         assert solve_params['n_s'] == 2
-        alpha = np.zeros((n_terms, 2, 2, 2))
+        alpha = np.zeros((n_terms + int(use_D) + int(use_Jperp), 2, 2, 2))
         for l, (term, coeff) in enumerate(h_int):
             bl0, bl1, u0, u0p, u1, u1p = self.indices_from_quartic_term(term, gf_struct)
 
@@ -309,6 +313,17 @@ def trivial_alpha(self, solve_params):
             else:
                 assert False, "I don't know this type of term"
 
+        if use_D:
+            # F.F. Assaad, T.C. Lang, Phys. Rev. B 76, 035116 (2007) -- Eq. (36)
+            alpha_s = lambda s: np.array([[ 0.5 + s*delta[0], 0.0              ],
+                                          [ 0.0             , 0.5 + s*delta[0] ]])
+            alpha[n_terms,...,0] = alpha_s(+1)
+            alpha[n_terms,...,1] = alpha_s(-1)
+
+        if use_Jperp:
+            alpha[n_terms + int(use_D),...,0] = 0.0
+            alpha[n_terms + int(use_D),...,1] = 0.0
+
         return alpha
 
     def solve(self, **solve_params):

test/python/dynamic/all.py

@@ -51,6 +51,8 @@
 
 # --------- Solve! ----------
 S.solve(h_int=h_int,
+        alpha_mode = "trivial",
+        n_s = 2,
         n_cycles = n_cyc,
         length_cycle = 50,
         n_warmup_cycles = 100,

test/python/dynamic/densdens.py

@@ -46,6 +46,8 @@
 
 # --------- Solve! ----------
 S.solve(h_int=h_int,
+        alpha_mode = "trivial",
+        n_s = 2,
         n_cycles = n_cyc,
         length_cycle = 50,
         n_warmup_cycles = 100,

test/python/dynamic/jperp.py

@@ -46,6 +46,8 @@
 
 # --------- Solve! ----------
 S.solve(h_int=h_int,
+        alpha_mode = "trivial",
+        n_s = 2,
         n_cycles = n_cyc,
         length_cycle = 50,
         n_warmup_cycles = 100,