diff --git a/c++/triqs_ctint/solver_core.cpp b/c++/triqs_ctint/solver_core.cpp index 32e0233c..53382538 100644 --- a/c++/triqs_ctint/solver_core.cpp +++ b/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{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{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; } } } diff --git a/c++/triqs_ctint/vertex_factories.cpp b/c++/triqs_ctint/vertex_factories.cpp index 4a9350e5..a15ad4b6 100644 --- a/c++/triqs_ctint/vertex_factories.cpp +++ b/c++/triqs_ctint/vertex_factories.cpp @@ -9,6 +9,10 @@ namespace triqs_ctint { std::vector vertex_factories; + int const n_terms = std::distance(params.h_int.begin(), params.h_int.end()); + bool const use_D = static_cast(D0_iw); + bool const use_Jperp = static_cast(Jperp_iw); + // ------------ Create Vertex Factory for Static Interactions -------------- { std::vector 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 indices; #ifdef INTERACTION_IS_COMPLEX @@ -80,14 +84,14 @@ 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); @@ -95,7 +99,7 @@ namespace triqs_ctint { } // ------------ Create Vertex Factory for Dynamic Spin-Spin Interactions -------------- - if (Jperp_iw) { + if (use_Jperp) { std::vector 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); diff --git a/python/triqs_ctint/solver.py b/python/triqs_ctint/solver.py index d9048cac..fe8c699a 100644 --- a/python/triqs_ctint/solver.py +++ b/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): diff --git a/test/python/dynamic/all.py b/test/python/dynamic/all.py index 33d457d4..aad60630 100644 --- a/test/python/dynamic/all.py +++ b/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, diff --git a/test/python/dynamic/densdens.py b/test/python/dynamic/densdens.py index 62bc1557..7ed91e8b 100644 --- a/test/python/dynamic/densdens.py +++ b/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, diff --git a/test/python/dynamic/jperp.py b/test/python/dynamic/jperp.py index 6197a7d4..71beafa6 100644 --- a/test/python/dynamic/jperp.py +++ b/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,