Skip to content

Commit be85e5d

Browse files
paco-ripaco-ri
committed
Move third-order comp'ns to examples dir
1 parent 84b7eb6 commit be85e5d

4 files changed

Lines changed: 605 additions & 0 deletions

File tree

‎CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -137,6 +137,7 @@ add_subdirectory(c++/${PROJECT_NAME})
137137
# Tests
138138
if(Build_Tests)
139139
add_subdirectory(test)
140+
add_subdirectory(examples)
140141
endif()
141142

142143
# Python

‎examples/CMakeLists.txt‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
# List of all examples
2+
set(all_examples
3+
dense_backbone_third_order.cpp
4+
block_sparse_backbone_third_order.cpp
5+
)
6+
7+
foreach(example ${all_examples})
8+
get_filename_component(example_name ${example} NAME_WE)
9+
get_filename_component(example_dir ${example} DIRECTORY)
10+
add_executable(${example_name} ${example})
11+
target_link_libraries(${example_name} ${PROJECT_NAME}::${PROJECT_NAME}_c ${PROJECT_NAME}_warnings)
12+
set_property(TARGET ${example_name} PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/${example_dir})
13+
# Run clang-tidy if found
14+
if(CLANG_TIDY_EXECUTABLE)
15+
set_target_properties(${example_name} PROPERTIES CXX_CLANG_TIDY "${CLANG_TIDY_EXECUTABLE}")
16+
endif()
17+
# Run cppcheck if found
18+
if(CPPCHECK_EXECUTABLE)
19+
add_custom_command(
20+
TARGET ${example_name}
21+
COMMAND ${CPPCHECK_EXECUTABLE}
22+
--enable=warning,style,performance,portability
23+
--std=c++20
24+
--template=gcc
25+
--verbose
26+
--force
27+
--quiet
28+
${CMAKE_CURRENT_SOURCE_DIR}/${example}
29+
)
30+
endif()
31+
endforeach()

‎examples/block_sparse_backbone_third_order.cpp‎

Lines changed: 318 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,318 @@
1+
#include <chrono>
2+
#include <nda/nda.hpp>
3+
#include <cppdlr/cppdlr.hpp>
4+
#include <triqs_xca/block_sparse.hpp>
5+
#include <triqs_xca/block_sparse_manual.hpp>
6+
#include <triqs_xca/block_sparse_backbone.hpp>
7+
8+
using namespace nda;
9+
using namespace cppdlr;
10+
11+
nda::array<dcomplex, 3> Hmat_to_Gtmat(nda::array<dcomplex, 2> Hmat, double beta, nda::array<double, 1> dlr_it_abs) {
12+
// Helper function for computing the non-interacting Green's function from the Hamiltonian, both in dense storage
13+
14+
int N = Hmat.extent(0);
15+
auto [H_loc_eval, H_loc_evec] = nda::linalg::eigenelements(Hmat);
16+
auto E0 = nda::min_element(H_loc_eval);
17+
H_loc_eval -= E0;
18+
auto tr_exp_minusbetaH = nda::sum(exp(-beta * H_loc_eval));
19+
auto eta_0 = nda::log(tr_exp_minusbetaH) / beta;
20+
H_loc_eval += eta_0;
21+
auto Gt_evals_t = nda::zeros<dcomplex>(N, N);
22+
int r = dlr_it_abs.extent(0);
23+
auto Gt_mat = nda::zeros<dcomplex>(r, N, N);
24+
auto Gbeta = nda::zeros<dcomplex>(N, N);
25+
for (int t = 0; t < r; t++) {
26+
for (int i = 0; i < N; i++) { Gt_evals_t(i, i) = -exp(-beta * dlr_it_abs(t) * H_loc_eval(i)); }
27+
Gt_mat(t, _, _) = matmul(H_loc_evec, matmul(Gt_evals_t, nda::transpose(H_loc_evec)));
28+
}
29+
return Gt_mat;
30+
}
31+
32+
std::tuple<nda::array<dcomplex, 3>, nda::array<dcomplex, 3>> discrete_bath_helper(double beta, double Lambda, double eps) {
33+
// Helper function for setting up the discrete bath model
34+
35+
auto dlr_rf = build_dlr_rf(Lambda, eps);
36+
auto itops = imtime_ops(Lambda, dlr_rf);
37+
auto const &dlr_it = itops.get_itnodes();
38+
auto dlr_it_abs = rel2abs(dlr_it);
39+
int r = itops.rank();
40+
41+
// hybridization parameters
42+
double s = 0.5;
43+
double t = 1.0;
44+
nda::array<double, 1> e{-2.3 * t, 2.3 * t};
45+
46+
// hybridization generation
47+
auto Jt = nda::array<dcomplex, 3>(r, 1, 1);
48+
auto Jt_refl = nda::array<dcomplex, 3>(r, 1, 1);
49+
for (int i = 0; i <= 1; i++) {
50+
for (int u = 0; u < r; u++) {
51+
Jt(u, 0, 0) += k_it(dlr_it(u), e(i), beta);
52+
Jt_refl(u, 0, 0) += k_it(-dlr_it(u), e(i), beta);
53+
}
54+
}
55+
56+
// orbital index order: do 0, do 1, up 0, up 1. same level <-> same parity index
57+
auto Deltat = nda::array<dcomplex, 3>(r, 4, 4);
58+
auto Deltat_refl = nda::array<dcomplex, 3>(r, 4, 4);
59+
60+
for (int i = 0; i < Deltat.extent(1); i++) {
61+
for (int j = 0; j < Deltat.extent(2); j++) {
62+
if (i == j) {
63+
Deltat(_, i, j) = Jt(_, 0, 0);
64+
Deltat_refl(_, i, j) = Jt_refl(_, 0, 0);
65+
} else if ((i == 0 && j == 1) || (i == 1 && j == 0) || (i == 2 && j == 3) || (i == 3 && j == 2)) {
66+
Deltat(_, i, j) = s * Jt(_, 0, 0);
67+
Deltat_refl(_, i, j) = s * Jt_refl(_, 0, 0);
68+
}
69+
}
70+
}
71+
Deltat = t * t * Deltat;
72+
Deltat_refl = t * t * Deltat_refl;
73+
74+
return std::make_tuple(Deltat, Deltat_refl);
75+
}
76+
77+
std::tuple<nda::array<dcomplex, 3>, nda::array<dcomplex, 3>, nda::array<dcomplex, 3>> two_band_dense_helper(double beta, double Lambda, double eps) {
78+
79+
auto dlr_rf = build_dlr_rf(Lambda, eps);
80+
auto itops = imtime_ops(Lambda, dlr_rf);
81+
auto const &dlr_it = itops.get_itnodes();
82+
auto dlr_it_abs = cppdlr::rel2abs(dlr_it);
83+
84+
// Hamiltonian in dense storage
85+
auto H_dense = nda::zeros<dcomplex>(16, 16);
86+
H_dense(0, 0) = -1;
87+
H_dense(1, 1) = -1;
88+
H_dense(2, 2) = -1;
89+
H_dense(3, 3) = -1;
90+
H_dense(4, 4) = -0.6;
91+
H_dense(5, 8) = 0.2;
92+
H_dense(6, 6) = -0.4;
93+
H_dense(6, 7) = 0.2;
94+
H_dense(7, 6) = 0.2;
95+
H_dense(7, 7) = -0.4;
96+
H_dense(8, 5) = 0.2;
97+
H_dense(9, 9) = -0.6;
98+
H_dense(11, 11) = 2;
99+
H_dense(12, 12) = 2;
100+
H_dense(13, 13) = 2;
101+
H_dense(14, 14) = 2;
102+
H_dense(15, 15) = 6;
103+
104+
// Green's function in dense storage
105+
auto Gt_dense = Hmat_to_Gtmat(H_dense, beta, dlr_it_abs);
106+
107+
// creation/annihilation operators in dense storage
108+
auto Fs_dense = nda::zeros<dcomplex>(4, 16, 16);
109+
// copied from a text dump of an h5 file output from atom_diag
110+
Fs_dense = {{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
111+
{0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
112+
{0, 0, 0, 0, 0, 1, 0, 0, 2.23711e-17, 0, 0, 0, 0, 0, 0, 0},
113+
{0, 0, 0, 0, 0, 0, 2.23711e-17, 1, 0, 0, 0, 0, 0, 0, 0, 0},
114+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
115+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.23711e-17, 0, 0, 0},
116+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
117+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.23711e-17, 0, 0, 0, 0},
118+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
119+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
120+
{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
121+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
122+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
123+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
124+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1},
125+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
126+
{{0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
127+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
128+
{0, 0, 0, 0, 0, 0, 1, 2.23711e-17, 0, 0, 0, 0, 0, 0, 0, 0},
129+
{0, 0, 0, 0, 0, 2.23711e-17, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
130+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
131+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0},
132+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.23711e-17, 0, 0, 0},
133+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0},
134+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.23711e-17, 0, 0, 0, 0},
135+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
136+
{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
137+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
138+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
139+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1},
140+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
141+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
142+
{{0, 0, 0, 0, 0, -1, 0, 0, -2.23711e-17, 0, 0, 0, 0, 0, 0, 0},
143+
{0, 0, 0, 0, 0, 0, -1, -2.23711e-17, 0, 0, 0, 0, 0, 0, 0, 0},
144+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
145+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0},
146+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
147+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.23711e-17, 0},
148+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.23711e-17, 0, 0},
149+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0},
150+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0},
151+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
152+
{0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
153+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
154+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1},
155+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
156+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
157+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
158+
{{0, 0, 0, 0, 0, 0, -2.23711e-17, -1, 0, 0, 0, 0, 0, 0, 0, 0},
159+
{0, 0, 0, 0, 0, -2.23711e-17, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0},
160+
{0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0},
161+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
162+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
163+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0},
164+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
165+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.23711e-17, 0},
166+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.23711e-17, 0, 0},
167+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
168+
{0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
169+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1},
170+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
171+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
172+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
173+
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}};
174+
auto F_dags_dense = nda::zeros<dcomplex>(4, 16, 16);
175+
for (int i = 0; i < 4; i++) { F_dags_dense(i, _, _) = nda::transpose(nda::conj(Fs_dense(i, _, _))); }
176+
177+
return std::make_tuple(Gt_dense, Fs_dense, F_dags_dense);
178+
}
179+
180+
std::tuple<BlockDiagOpFun, BlockOpSymQuartet, nda::vector<int>> two_band_helper(double beta, double Lambda, double eps,
181+
nda::array_const_view<dcomplex, 3> hyb_coeffs,
182+
nda::array_const_view<dcomplex, 3> hyb_refl_coeffs) {
183+
auto dlr_rf = build_dlr_rf(Lambda, eps);
184+
auto itops = imtime_ops(Lambda, dlr_rf);
185+
auto const &dlr_it = itops.get_itnodes();
186+
auto dlr_it_abs = cppdlr::rel2abs(dlr_it);
187+
188+
// get Hamiltonian, creation/annihilation operators in block-sparse storage
189+
int num_blocks = 5; // number of blocks of Hamiltonian
190+
191+
// Hamiltonian
192+
std::vector<nda::array<double, 2>> H_blocks(num_blocks); // Hamiltonian in sparse storage
193+
H_blocks[0] = nda::make_regular(-1 * nda::eye<double>(4));
194+
H_blocks[1] = {{-0.6, 0, 0, 0, 0, 0}, {0, 8.27955e-19, 0, 0, 0.2, 0}, {0, 0, -0.4, 0.2, 0, 0},
195+
{0, 0, 0.2, -0.4, 0, 0}, {0, 0.2, 0, 0, 8.27955e-19, 0}, {0, 0, 0, 0, 0, -0.6}};
196+
H_blocks[2] = {{0}};
197+
H_blocks[3] = nda::make_regular(2 * nda::eye<double>(4));
198+
H_blocks[4] = {{6}};
199+
nda::vector<int> H_block_inds = {0, 0, -1, 0, 0};
200+
201+
// Green's function
202+
auto Gt = nonint_gf_BDOF(H_blocks, H_block_inds, beta, dlr_it_abs);
203+
204+
// creation/annihilation operators
205+
nda::vector<int> ann_conn = {2, 0, -1, 1, 3}; // block column indices of F operators
206+
nda::vector<int> cre_conn = {1, 3, 0, 4, -1}; // block column indices of F^dag operators
207+
std::vector<nda::array<dcomplex, 3>> F_blocks(num_blocks), Fdag_blocks(num_blocks);
208+
209+
F_blocks[0] = {{{1, 0, 0, 0}}, {{0, 1, 0, 0}}, {{0, 0, 1, 0}}, {{0, 0, 0, 1}}};
210+
F_blocks[1] = {{{0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0}},
211+
{{-1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 1, 0}},
212+
{{0, -1, 0, 0, 0, 0}, {0, 0, -1, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1}},
213+
{{0, 0, 0, -1, 0, 0}, {0, 0, 0, 0, -1, 0}, {0, 0, 0, 0, 0, -1}, {0, 0, 0, 0, 0, 0}}};
214+
F_blocks[2] = {{{0}}};
215+
F_blocks[3] = {{{0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}},
216+
{{0, 0, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}},
217+
{{1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}, {0, 0, 0, 0}},
218+
{{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}};
219+
F_blocks[4] = {{{0}, {0}, {0}, {1}}, {{0}, {0}, {-1}, {0}}, {{0}, {1}, {0}, {0}}, {{-1}, {0}, {0}, {0}}};
220+
221+
Fdag_blocks[0] = {{{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}},
222+
{{-1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}},
223+
{{0, 0, 0, 0}, {-1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}},
224+
{{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}}};
225+
Fdag_blocks[1] = {{{0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, 0}},
226+
{{0, -1, 0, 0, 0, 0}, {0, 0, 0, -1, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1}},
227+
{{1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, -1, 0, 0}, {0, 0, 0, 0, -1, 0}},
228+
{{0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0}}};
229+
Fdag_blocks[2] = {{{1}, {0}, {0}, {0}}, {{0}, {1}, {0}, {0}}, {{0}, {0}, {1}, {0}}, {{0}, {0}, {0}, {1}}};
230+
Fdag_blocks[3] = {{{0, 0, 0, 1}}, {{0, 0, -1, 0}}, {{0, 1, 0, 0}}, {{-1, 0, 0, 0}}};
231+
Fdag_blocks[4] = {{{0}}};
232+
233+
BlockOpSymSet F_BOSS(ann_conn, F_blocks), Fdag_BOSS(cre_conn, Fdag_blocks);
234+
std::vector<BlockOpSymSet> F_sym_vec{F_BOSS}, F_dag_sym_vec{Fdag_BOSS};
235+
nda::vector<long> sym_set_labels(4);
236+
sym_set_labels = 0; // all operators belong to the same symmetry set
237+
BlockOpSymQuartet Fq(F_sym_vec, F_dag_sym_vec, hyb_coeffs, hyb_refl_coeffs, sym_set_labels);
238+
239+
return std::make_tuple(Gt, Fq, sym_set_labels);
240+
}
241+
242+
int main() {
243+
nda::array<int, 3> topologies = {{{0, 2}, {1, 4}, {3, 5}}, {{0, 3}, {1, 5}, {2, 4}}, {{0, 4}, {1, 3}, {2, 5}}, {{0, 3}, {1, 4}, {2, 5}}};
244+
nda::vector<int> topo_sign{1, 1, 1, -1}; // topo_sign(i) = (-1)^{# of line crossings in topology i}
245+
246+
nda::array<int, 2> topology = {{0, 2}, {1, 3}};
247+
int n = 4, N = 16;
248+
double beta = 2.0;
249+
double Lambda = 10.0 * beta;
250+
double eps = 1.0e-4;
251+
252+
// generate hybridization, noninteracting Green's function, creation/annihilation operators
253+
auto [Deltat, Deltat_refl] = discrete_bath_helper(beta, Lambda, eps);
254+
auto [Gt_dense, Fs_dense, F_dags_dense] = two_band_dense_helper(beta, Lambda, eps);
255+
256+
// DLR generation
257+
auto dlr_rf = build_dlr_rf(Lambda, eps);
258+
auto itops = imtime_ops(Lambda, dlr_rf);
259+
auto const &dlr_it = itops.get_itnodes();
260+
auto dlr_it_abs = cppdlr::rel2abs(dlr_it);
261+
int r = itops.rank();
262+
263+
// hybridization and DenseFSet
264+
auto hyb_coeffs = itops.vals2coefs(Deltat); // hybridization DLR coeffs
265+
auto hyb_refl = Deltat;
266+
auto hyb_refl_coeffs = hyb_coeffs;
267+
auto [Gt, Fq, sym_set_labels] = two_band_helper(beta, Lambda, eps, hyb_coeffs, hyb_refl_coeffs);
268+
269+
auto D = DiagramBlockSparseEvaluator(beta, itops, Deltat, hyb_refl, dlr_rf, Gt, Fq); // create DiagramEvaluator object
270+
271+
auto Deltadlr = itops.vals2coefs(Deltat); //obtain dlr coefficient of Delta(t)
272+
nda::vector<double> dlr_rf_reflect = -dlr_rf;
273+
nda::array<dcomplex, 3> Deltadlr_reflect = Deltadlr * 1.0;
274+
for (int i = 0; i < Deltadlr.shape(0); ++i) Deltadlr_reflect(i, _, _) = transpose(Deltadlr(i, _, _));
275+
auto Delta_decomp = hyb_decomp(Deltadlr, dlr_rf, eps);
276+
auto Delta_decomp_reflect = hyb_decomp(Deltadlr_reflect, dlr_rf_reflect, eps);
277+
hyb_F Delta_F(N, r, n);
278+
hyb_F Delta_F_reflect(N, r, n);
279+
Delta_F.update_inplace(Delta_decomp, dlr_it, Fs_dense, F_dags_dense);
280+
Delta_F_reflect.update_inplace(Delta_decomp_reflect, dlr_it, F_dags_dense, Fs_dense);
281+
auto fb = nda::vector<int64_t>(3);
282+
fb(0) = 0;
283+
284+
BlockDiagOpFun third_order_result(r, Gt.get_block_sizes());
285+
for (int i = 0; i < 4; ++i) {
286+
std::cout << "Evaluating topology " << i << std::endl;
287+
288+
// Compute third-order contribution using DiagramEvaluator
289+
auto B = Backbone(topologies(i, _, _), n);
290+
D.eval_diagram_block_sparse(B);
291+
third_order_result = D.Sigma;
292+
D.reset(); // reset the DiagramEvaluator for the next topology
293+
294+
// Compute third-order contribution using Zhen's code
295+
nda::array<dcomplex, 3> TCA_Zhen(r, N, N);
296+
TCA_Zhen = 0;
297+
fb(1) = 0;
298+
fb(2) = 0;
299+
TCA_Zhen += Sigma_Diagram_calc(Delta_F, Delta_F_reflect, topologies(i, _, _), Deltat, Deltat_refl, Gt_dense, itops, beta, Fs_dense,
300+
F_dags_dense, fb, true);
301+
fb(1) = 1;
302+
TCA_Zhen += Sigma_Diagram_calc(Delta_F, Delta_F_reflect, topologies(i, _, _), Deltat, Deltat_refl, Gt_dense, itops, beta, Fs_dense,
303+
F_dags_dense, fb, true);
304+
fb(2) = 1;
305+
TCA_Zhen += Sigma_Diagram_calc(Delta_F, Delta_F_reflect, topologies(i, _, _), Deltat, Deltat_refl, Gt_dense, itops, beta, Fs_dense,
306+
F_dags_dense, fb, true);
307+
fb(1) = 0;
308+
TCA_Zhen += Sigma_Diagram_calc(Delta_F, Delta_F_reflect, topologies(i, _, _), Deltat, Deltat_refl, Gt_dense, itops, beta, Fs_dense,
309+
F_dags_dense, fb, true);
310+
311+
std::cout << "max error in block 0: " << nda::max_element(nda::abs(TCA_Zhen(_, range(0, 4), range(0, 4)) - third_order_result.get_block(0))) << std::endl;
312+
std::cout << "max error in block 1: " << nda::max_element(nda::abs(TCA_Zhen(_, range(4, 10), range(4, 10)) - third_order_result.get_block(1))) << std::endl;
313+
std::cout << "max error in block 2: " << nda::max_element(nda::abs(TCA_Zhen(_, range(10, 11), range(10, 11)) - third_order_result.get_block(2))) << std::endl;
314+
std::cout << "max error in block 3: " << nda::max_element(nda::abs(TCA_Zhen(_, range(11, 15), range(11, 15)) - third_order_result.get_block(3))) << std::endl;
315+
std::cout << "max error in block 4: " << nda::max_element(nda::abs(TCA_Zhen(_, range(15, 16), range(15, 16)) - third_order_result.get_block(4))) << std::endl;
316+
std::cout << std::endl;
317+
}
318+
}

0 commit comments

Comments
 (0)