TRIQS/triqs_ctint 4.0.0
A TRIQS application
Loading...
Searching...
No Matches
post_process.cpp
1// Copyright (c) 2017--present, The Simons Foundation
2// This file is part of TRIQS/ctint and is licensed under the terms of GPLv3 or later.
3// SPDX-License-Identifier: GPL-3.0-or-later
4// See LICENSE in the root of this distribution for details.
5
6#include "post_process.hpp"
7#include <triqs/gfs.hpp>
8#include <triqs/mesh.hpp>
9
10namespace triqs_ctint {
11
12 chi4_iw_t G2_conn_from_M4(chi4_iw_t::const_view_type M4_iw, g_reg_iw_t::const_view_type M_iw, g_reg_iw_t::const_view_type G0_iw) {
13
14 chi4_iw_t G2_conn_iw = M4_iw; // FIXME Product Ranges with += Lazy Expressions
15
16 double beta = M_iw[0].mesh().beta();
17 int n_blocks = M_iw.size();
18
19 // Calculate connected part of M4
20 chi4_iw_t M4_iw_conn = M4_iw;
21
22 for (int bl1 : range(n_blocks))
23 for (int bl2 : range(n_blocks))
24 M4_iw_conn(bl1, bl2)(iw1_, iw2_, iw3_)(i_, j_, k_, l_) << M4_iw(bl1, bl2)(iw1_, iw2_, iw3_)(i_, j_, k_, l_)
25 - beta * kronecker(iw1_, iw2_) * M_iw[bl1](iw1_)(j_, i_) * M_iw[bl2](iw3_)(l_, k_)
26 + beta * kronecker(bl1, bl2) * kronecker(iw2_, iw3_) * M_iw[bl1](iw1_)(l_, i_) * M_iw[bl2](iw3_)(j_, k_);
27
28 // Calculate disconnected part of the two-particle Green function
29 G2_conn_iw() = 0.;
30
31 for (int bl1 : range(n_blocks))
32 for (int bl2 : range(n_blocks)) {
33
34 int bl1_size = M4_iw(bl1, bl2).target_shape()[0];
35 int bl2_size = M4_iw(bl1, bl2).target_shape()[2];
36
37 for (int m : range(bl1_size))
38 for (int n : range(bl1_size))
39 for (int o : range(bl2_size))
40 for (int p : range(bl2_size))
41 G2_conn_iw(bl1, bl2)(iw1_, iw2_, iw3_)(i_, j_, k_, l_) << G2_conn_iw(bl1, bl2)(iw1_, iw2_, iw3_)(i_, j_, k_, l_)
42 + G0_iw[bl1](iw2_)(j_, n) * G0_iw[bl2](iw1_ - iw2_ + iw3_)(l_, p) * M4_iw_conn(bl1, bl2)(iw1_, iw2_, iw3_)(m, n, o, p)
43 * G0_iw[bl1](iw1_)(m, i_) * G0_iw[bl2](iw3_)(o, k_);
44 }
45
46 return G2_conn_iw;
47 }
48
49 chi4_iw_t G2pp_conn_from_M4pp(chi4_iw_t::const_view_type M4pp_iw, g_reg_iw_t::const_view_type M_iw, g_reg_iw_t::const_view_type G0_iw) {
50
51 chi4_iw_t G2pp_conn_iw = M4pp_iw; // FIXME Product Ranges with += Lazy Expressions
52
53 double beta = M_iw[0].mesh().beta();
54 int n_blocks = M_iw.size();
55
56 // Calculate connected part of M4
57 chi4_iw_t M4pp_iw_conn = M4pp_iw;
58
59 for (int bl1 : range(n_blocks))
60 for (int bl2 : range(n_blocks))
61 M4pp_iw_conn(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << M4pp_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
62 - beta * kronecker(iw_, iW_ - iwp_) * M_iw[bl1](iw_)(j_, i_) * M_iw[bl2](iW_ - iw_)(l_, k_)
63 + beta * kronecker(bl1, bl2) * kronecker(iW_ - iwp_, iW_ - iw_) * M_iw[bl1](iw_)(l_, i_) * M_iw[bl2](iW_ - iw_)(j_, k_);
64
65 // Calculate disconnected part of the two-particle Green function
66 G2pp_conn_iw() = 0.;
67
68 for (int bl1 : range(n_blocks))
69 for (int bl2 : range(n_blocks)) {
70
71 int bl1_size = M4pp_iw(bl1, bl2).target_shape()[0];
72 int bl2_size = M4pp_iw(bl1, bl2).target_shape()[2];
73
74 for (int m : range(bl1_size))
75 for (int n : range(bl1_size))
76 for (int o : range(bl2_size))
77 for (int p : range(bl2_size))
78 G2pp_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << G2pp_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
79 + G0_iw[bl1](iW_ - iwp_)(j_, n) * G0_iw[bl2](iwp_)(l_, p) * M4pp_iw_conn(bl1, bl2)(iW_, iw_, iwp_)(m, n, o, p)
80 * G0_iw[bl1](iw_)(m, i_) * G0_iw[bl2](iW_ - iw_)(o, k_);
81 }
82
83 return G2pp_conn_iw;
84 }
85
86 chi4_iw_t G2ph_conn_from_M4ph(chi4_iw_t::const_view_type M4ph_iw, g_reg_iw_t::const_view_type M_iw, g_reg_iw_t::const_view_type G0_iw) {
87
88 chi4_iw_t G2ph_conn_iw = M4ph_iw; // FIXME Product Ranges with += Lazy Expressions
89
90 double beta = M_iw[0].mesh().beta();
91 int n_blocks = M_iw.size();
92
93 // Calculate connected part of M4
94 chi4_iw_t M4ph_iw_conn = M4ph_iw;
95
96 for (int bl1 : range(n_blocks))
97 for (int bl2 : range(n_blocks))
98 M4ph_iw_conn(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << M4ph_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
99 - beta * kronecker(iw_, iW_ + iw_) * M_iw[bl1](iw_)(j_, i_) * M_iw[bl2](iW_ + iwp_)(l_, k_)
100 + beta * kronecker(bl1, bl2) * kronecker(iW_ + iw_, iW_ + iwp_) * M_iw[bl1](iw_)(l_, i_) * M_iw[bl2](iW_ + iwp_)(j_, k_);
101
102 // Calculate disconnected part of the two-particle Green function
103 G2ph_conn_iw() = 0.;
104
105 for (int bl1 : range(n_blocks))
106 for (int bl2 : range(n_blocks)) {
107
108 int bl1_size = M4ph_iw(bl1, bl2).target_shape()[0];
109 int bl2_size = M4ph_iw(bl1, bl2).target_shape()[2];
110
111 for (int m : range(bl1_size))
112 for (int n : range(bl1_size))
113 for (int o : range(bl2_size))
114 for (int p : range(bl2_size))
115 G2ph_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << G2ph_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
116 + G0_iw[bl1](iW_ + iw_)(j_, n) * G0_iw[bl2](iwp_)(l_, p) * M4ph_iw_conn(bl1, bl2)(iW_, iw_, iwp_)(m, n, o, p)
117 * G0_iw[bl1](iw_)(m, i_) * G0_iw[bl2](iW_ + iwp_)(o, k_);
118 }
119
120 return G2ph_conn_iw;
121 }
122
123 chi4_iw_t F_from_G2c(chi4_iw_t::const_view_type G2_conn_iw, g_reg_iw_t::const_view_type G_iw) {
124
125 int n_blocks = G_iw.size();
126
127 // Temporary quantities
128 g_reg_iw_t Ginv = inverse(G_iw);
129
130 // Calculate vertex function F
131 chi4_iw_t F_iw = G2_conn_iw; // FIXME Product Ranges with += Lazy Expressions
132 F_iw() = 0;
133
134 for (int bl1 : range(n_blocks))
135 for (int bl2 : range(n_blocks)) {
136
137 int bl1_size = G2_conn_iw(bl1, bl2).target_shape()[0];
138 int bl2_size = G2_conn_iw(bl1, bl2).target_shape()[2];
139
140 for (int m : range(bl1_size))
141 for (int n : range(bl1_size))
142 for (int o : range(bl2_size))
143 for (int p : range(bl2_size))
144 F_iw(bl1, bl2)(iw1_, iw2_, iw3_)(i_, j_, k_, l_) << F_iw(bl1, bl2)(iw1_, iw2_, iw3_)(i_, j_, k_, l_)
145 + Ginv[bl1](iw2_)(j_, n) * Ginv[bl2](iw1_ - iw2_ + iw3_)(l_, p) * G2_conn_iw(bl1, bl2)(iw1_, iw2_, iw3_)(m, n, o, p)
146 * Ginv[bl1](iw1_)(m, i_) * Ginv[bl2](iw3_)(o, k_);
147 }
148
149 return F_iw;
150 }
151
152 chi4_iw_t Fpp_from_G2pp_conn(chi4_iw_t::const_view_type G2pp_conn_iw, g_reg_iw_t::const_view_type G_iw) {
153
154 int n_blocks = G_iw.size();
155
156 // Temporary quantities
157 g_reg_iw_t Ginv = inverse(G_iw);
158
159 // Calculate vertex function F
160 chi4_iw_t Fpp_iw = G2pp_conn_iw; // FIXME Product Ranges with += Lazy Expressions
161 Fpp_iw() = 0;
162
163 for (int bl1 : range(n_blocks))
164 for (int bl2 : range(n_blocks)) {
165
166 int bl1_size = G2pp_conn_iw(bl1, bl2).target_shape()[0];
167 int bl2_size = G2pp_conn_iw(bl1, bl2).target_shape()[2];
168
169 for (int m : range(bl1_size))
170 for (int n : range(bl1_size))
171 for (int o : range(bl2_size))
172 for (int p : range(bl2_size))
173 Fpp_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << Fpp_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
174 + Ginv[bl1](iW_ - iwp_)(j_, n) * Ginv[bl2](iwp_)(l_, p) * G2pp_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(m, n, o, p) * Ginv[bl1](iw_)(m, i_)
175 * Ginv[bl2](iW_ - iw_)(o, k_);
176 }
177
178 return Fpp_iw;
179 }
180
181 chi4_iw_t Fph_from_G2ph_conn(chi4_iw_t::const_view_type G2ph_conn_iw, g_reg_iw_t::const_view_type G_iw) {
182
183 int n_blocks = G_iw.size();
184
185 // Temporary quantities
186 g_reg_iw_t Ginv = inverse(G_iw);
187
188 // Calculate vertex function F
189 chi4_iw_t Fph_iw = G2ph_conn_iw; // FIXME Product Ranges with += Lazy Expressions
190 Fph_iw() = 0;
191
192 for (int bl1 : range(n_blocks))
193 for (int bl2 : range(n_blocks)) {
194
195 int bl1_size = G2ph_conn_iw(bl1, bl2).target_shape()[0];
196 int bl2_size = G2ph_conn_iw(bl1, bl2).target_shape()[2];
197
198 for (int m : range(bl1_size))
199 for (int n : range(bl1_size))
200 for (int o : range(bl2_size))
201 for (int p : range(bl2_size))
202 Fph_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << Fph_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
203 + Ginv[bl1](iW_ + iw_)(j_, n) * Ginv[bl2](iwp_)(l_, p) * G2ph_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(m, n, o, p) * Ginv[bl1](iw_)(m, i_)
204 * Ginv[bl2](iW_ + iwp_)(o, k_);
205 }
206
207 return Fph_iw;
208 }
209
210 chi4_iw_t G2_from_G2c(chi4_iw_t::const_view_type G2_conn_iw, g_reg_iw_t::const_view_type G_iw) {
211
212 int n_blocks = G_iw.size();
213 double beta = G_iw[0].mesh().beta();
214
215 // Calculate G2_iw from G2_conn_iw and G_iw
216 chi4_iw_t G2_iw = G2_conn_iw;
217
218 for (int bl1 : range(n_blocks))
219 for (int bl2 : range(n_blocks))
220 G2_iw(bl1, bl2)(iw1_, iw2_, iw3_)(i_, j_, k_, l_) << G2_conn_iw(bl1, bl2)(iw1_, iw2_, iw3_)(i_, j_, k_, l_)
221 + beta * kronecker(iw1_, iw2_) * G_iw[bl1](iw1_)(j_, i_) * G_iw[bl2](iw3_)(l_, k_)
222 - beta * kronecker(bl1, bl2) * kronecker(iw2_, iw3_) * G_iw[bl1](iw1_)(l_, i_) * G_iw[bl2](iw3_)(j_, k_);
223
224 return G2_iw;
225 }
226
227 chi4_iw_t G2pp_from_G2pp_conn(chi4_iw_t::const_view_type G2pp_conn_iw, g_reg_iw_t::const_view_type G_iw) {
228
229 int n_blocks = G_iw.size();
230 double beta = G_iw[0].mesh().beta();
231
232 // Calculate G2_iw from G2_conn_iw and G_iw
233 chi4_iw_t G2pp_iw = G2pp_conn_iw;
234
235 for (int bl1 : range(n_blocks))
236 for (int bl2 : range(n_blocks))
237 G2pp_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << G2pp_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
238 + beta * kronecker(iw_, iW_ - iwp_) * G_iw[bl1](iw_)(j_, i_) * G_iw[bl2](iW_ - iw_)(l_, k_)
239 - beta * kronecker(bl1, bl2) * kronecker(iW_ - iwp_, iW_ - iw_) * G_iw[bl1](iw_)(l_, i_) * G_iw[bl2](iW_ - iw_)(j_, k_);
240
241 return G2pp_iw;
242 }
243
244 chi4_iw_t G2ph_from_G2ph_conn(chi4_iw_t::const_view_type G2ph_conn_iw, g_reg_iw_t::const_view_type G_iw) {
245
246 int n_blocks = G_iw.size();
247 double beta = G_iw[0].mesh().beta();
248
249 // Calculate G2_iw from G2_conn_iw and G_iw
250 chi4_iw_t G2ph_iw = G2ph_conn_iw;
251
252 for (int bl1 : range(n_blocks))
253 for (int bl2 : range(n_blocks))
254 G2ph_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << G2ph_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
255 + beta * kronecker(iw_, iW_ + iw_) * G_iw[bl1](iw_)(j_, i_) * G_iw[bl2](iW_ + iwp_)(l_, k_)
256 - beta * kronecker(bl1, bl2) * kronecker(iW_ + iw_, iW_ + iwp_) * G_iw[bl1](iw_)(l_, i_) * G_iw[bl2](iW_ + iwp_)(j_, k_);
257
258 return G2ph_iw;
259 }
260
261 chi4_iw_t chi_tilde_ph_from_G2ph_conn(chi4_iw_t::const_view_type G2ph_conn_iw, g_reg_iw_cv_t G_iw) {
262
263 int n_blocks = G_iw.size();
264 double beta = G_iw[0].mesh().beta();
265
266 // Calculate chi_tilde_pha from G2ph_conn_iw and G_iw
267 chi4_iw_t chi_tilde_ph = G2ph_conn_iw;
268
269 // Calculate generalized susceptibility in the ph channel from G2_conn_iw and G_iw
270 for (int bl1 : range(n_blocks))
271 for (int bl2 : range(n_blocks))
272 chi_tilde_ph(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_) << G2ph_conn_iw(bl1, bl2)(iW_, iw_, iwp_)(i_, j_, k_, l_)
273 - beta * kronecker(bl1, bl2) * kronecker(iw_, iwp_) * G_iw[bl1](iw_)(l_, i_) * G_iw[bl2](iwp_ + iW_)(j_, k_);
274
275 return chi_tilde_ph;
276 }
277
278} // namespace triqs_ctint