@@ -32,9 +32,9 @@ namespace triqs_ctseg::measures {
32
32
for (auto const &[bl1_name, bl1_size] : wdata.gf_struct )
33
33
for (auto const &[bl2_name, bl2_size] : wdata.gf_struct )
34
34
green_v.emplace_back (gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4 >>(
35
- {{beta, Fermion, p.n_iw_chi4_f , imfreq::option::all_frequencies},
36
- {beta, Fermion, p.n_iw_chi4_f , imfreq::option::all_frequencies},
37
- {beta, Boson, p.n_iw_chi4_b , imfreq::option::positive_frequencies_only}},
35
+ {{beta, Fermion, p.n_w_f_vertex , imfreq::option::all_frequencies},
36
+ {beta, Fermion, p.n_w_f_vertex , imfreq::option::all_frequencies},
37
+ {beta, Boson, p.n_w_b_vertex , imfreq::option::positive_frequencies_only}},
38
38
make_shape (bl1_size, bl1_size, bl2_size, bl2_size)));
39
39
return green_v;
40
40
};
@@ -48,18 +48,15 @@ namespace triqs_ctseg::measures {
48
48
49
49
g3w = make_block_gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4 >>(bosonic_block_names, g3w_vec ());
50
50
f3w = make_block_gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4 >>(bosonic_block_names, g3w_vec ());
51
-
52
- LOG ( " \n ====================== COMPUTING M-MATRIX ====================== \n " ) ;
51
+ g3w () = 0 ;
52
+ f3w () = 0 ;
53
53
54
54
n_w_fermionic = std::get<0 >(g3w[0 ].mesh ().components ()).last_index () + 1 ;
55
55
n_w_bosonic = std::get<2 >(g3w[0 ].mesh ().components ()).last_index () + 1 ;
56
56
57
57
w_ini = (2 * (-n_w_fermionic) + 1 ) * M_PI / beta;
58
58
w_inc = 2 * M_PI / beta;
59
59
60
- Mw_vector = compute_Mw (false );
61
- if (measure_f3w) nMw_vector = compute_Mw (true );
62
-
63
60
}
64
61
65
62
// -------------------------------------
@@ -89,6 +86,11 @@ namespace triqs_ctseg::measures {
89
86
90
87
Z += s;
91
88
89
+ auto const mesh_fermionic = std::get<0 >(g3w[0 ].mesh ());
90
+ auto const mesh_bosonic = std::get<2 >(g3w[0 ].mesh ());
91
+
92
+ Mw_vector = compute_Mw (false );
93
+
92
94
if (measure_g3w) {
93
95
for (long bl = 0 ; bl < g3w.size (); bl++) { // bl : 'upup', 'updn', ...
94
96
long b1 = bl / wdata.gf_struct .size ();
@@ -102,9 +104,13 @@ namespace triqs_ctseg::measures {
102
104
for (int m = 0 ; m < n_w_bosonic; m++) {
103
105
int n2 = n1 + m;
104
106
int n3 = n4 + m;
105
- g3w[bl](n1, n4, m)(a, b, c, d) = s * Mw (b1, a, b, n1, n2) * Mw (b2, c, d, n3, n4);
107
+ // This structure is ugly. Need someone who familiar with TRIQS to prune this part.
108
+ auto freq_1 = mesh_fermionic[n1 + n_w_fermionic];
109
+ auto freq_2 = mesh_fermionic[n4 + n_w_fermionic];
110
+ auto freq_3 = mesh_bosonic[m];
111
+ g3w[bl][freq_1, freq_2, freq_3](a, b, c, d) += s * Mw (b1, a, b, n1, n2) * Mw (b2, c, d, n3, n4);
106
112
if (b1 == b2)
107
- g3w[bl](n1, n4, m) (a, b, c, d) -= s * Mw (b1, a, d, n1, n4) * Mw (b2, c, b, n3, n2);
113
+ g3w[bl][freq_1, freq_2, freq_3] (a, b, c, d) -= s * Mw (b1, a, d, n1, n4) * Mw (b2, c, b, n3, n2);
108
114
} // m
109
115
} // n4
110
116
} // n1
@@ -116,6 +122,7 @@ namespace triqs_ctseg::measures {
116
122
} // measure_g3w
117
123
118
124
if (measure_f3w) {
125
+ nMw_vector = compute_Mw (true );
119
126
for (long bl = 0 ; bl < f3w.size (); bl++) { // bl : 'upup', 'updn', ...
120
127
long b1 = bl / wdata.gf_struct .size ();
121
128
long b2 = bl % wdata.gf_struct .size ();
@@ -128,9 +135,13 @@ namespace triqs_ctseg::measures {
128
135
for (int m = 0 ; m < n_w_bosonic; m++) {
129
136
int n2 = n1 + m;
130
137
int n3 = n4 + m;
131
- f3w[bl](n1, n4, m)(a, b, c, d) = s * nMw (b1, a, b, n1, n2) * Mw (b2, c, d, n3, n4);
138
+ // Please prune this
139
+ auto freq_1 = mesh_fermionic[n1 + n_w_fermionic];
140
+ auto freq_2 = mesh_fermionic[n4 + n_w_fermionic];
141
+ auto freq_3 = mesh_bosonic[m];
142
+ f3w[bl][freq_1, freq_2, freq_3](a, b, c, d) += s * nMw (b1, a, b, n1, n2) * Mw (b2, c, d, n3, n4);
132
143
if (b1 == b2)
133
- f3w[bl](n1, n4, m) (a, b, c, d) -= s * nMw (b1, a, d, n1, n4) * Mw (b2, c, b, n3, n2);
144
+ f3w[bl][freq_1, freq_2, freq_3] (a, b, c, d) -= s * nMw (b1, a, d, n1, n4) * Mw (b2, c, b, n3, n2);
134
145
} // m
135
146
} // n4
136
147
} // n1
@@ -196,8 +207,10 @@ namespace triqs_ctseg::measures {
196
207
auto const &[bl_name, bl_size] = bl_pair;
197
208
result[bl].resize (make_shape (bl_size, bl_size, n_w_aux, n_w_aux));
198
209
result[bl]() = 0 ;
210
+ }
199
211
200
- long N = wdata.dets [bl].size ();
212
+ for (auto const &[bl, det] : itertools::enumerate (wdata.dets )) {
213
+ long N = det.size ();
201
214
y_exp_ini.resize (N);
202
215
y_exp_inc.resize (N);
203
216
x_exp_ini.resize (N);
@@ -206,8 +219,8 @@ namespace triqs_ctseg::measures {
206
219
x_inner_index.resize (N);
207
220
208
221
for (long id : range (N)) {
209
- auto y = wdata. dets [bl] .get_y (id);
210
- auto x = wdata. dets [bl] .get_x (id);
222
+ auto y = det .get_y (id);
223
+ auto x = det .get_x (id);
211
224
y_exp_ini (id) = std::exp (dcomplex (0 , w_ini * double (std::get<0 >(y))));
212
225
y_exp_inc (id) = std::exp (dcomplex (0 , w_inc * double (std::get<0 >(y))));
213
226
x_exp_ini (id) = std::exp (dcomplex (0 , -w_ini * double (std::get<0 >(x))));
@@ -217,15 +230,15 @@ namespace triqs_ctseg::measures {
217
230
}
218
231
219
232
for (long id_y : range (N)) {
220
- auto y = wdata. dets [bl] .get_y (id_y);
233
+ auto y = det .get_y (id_y);
221
234
int yj = y_inner_index (id_y);
222
235
double f_fact = is_nMw ? fprefactor (bl, y) : 1.0 ;
223
236
224
237
for (long id_x : range (N)) {
225
238
int xi = x_inner_index (id_x);
226
239
dcomplex y_exp = y_exp_ini (id_y);
227
240
dcomplex x_exp = x_exp_ini (id_x);
228
- auto Minv = wdata. dets [bl] .inverse_matrix (id_y, id_x);
241
+ auto Minv = det .inverse_matrix (id_y, id_x);
229
242
230
243
for (int n_1 : range (n_w_aux)) {
231
244
for (int n_2 : range (n_w_aux)) {
0 commit comments