TRIQS/triqs_ctseg unstable
A TRIQS application
Loading...
Searching...
No Matches
four_point.cpp
1// Copyright (c) 2024--present, The Simons Foundation
2// Copyright (c) 2024--present, Max Planck Institute for Polymer Research, Mainz, Germany
3// This file is part of TRIQS/ctseg and is licensed under the terms of GPLv3 or later.
4// SPDX-License-Identifier: GPL-3.0-or-later
5// See LICENSE in the root of this distribution for details.
6
7#include "./four_point.hpp"
8#include "../logs.hpp"
9#include <triqs/utility/exceptions.hpp>
10
11namespace triqs_ctseg::measures {
12
13 four_point::four_point(params_t const &p, work_data_t const &wdata, configuration_t const &config, results_t &results)
14 : wdata{wdata}, config{config}, results{results} {
15
16 beta = p.beta;
17 n_w_bosonic = p.n_w_b_vertex;
18 n_w_fermionic = p.n_w_f_vertex;
19 measure_g2w = p.measure_g2w;
20 measure_g3w = p.measure_g3w;
21 mesh_bosonic = triqs::mesh::imfreq(beta, Boson , n_w_bosonic );
22 mesh_fermionic = triqs::mesh::imfreq(beta, Fermion, n_w_fermionic);
23
24 if (measure_g3w) {
25 if (p.length_cycle_g3 <= 0 or p.length_cycle_g3 % p.length_cycle != 0)
26 TRIQS_RUNTIME_ERROR << "length_cycle_g3 = " << p.length_cycle_g3
27 << " must be a positive multiple of length_cycle = " << p.length_cycle;
28 g3w_interval = p.length_cycle_g3 / p.length_cycle;
29 }
30
31 if (measure_g2w) g2w.resize(wdata.gf_struct.size());
32 if (measure_g3w) g3w.resize(wdata.gf_struct.size());
33 for (auto const &[bl1_idx, bl1] : itertools::enumerate(wdata.gf_struct)) {
34 auto &[bl1_name, bl1_size] = bl1;
35 block_names.push_back(bl1_name);
36 if (measure_g2w) g2w[bl1_idx].resize(wdata.gf_struct.size());
37 if (measure_g3w) g3w[bl1_idx].resize(wdata.gf_struct.size());
38 for (auto const &[bl2_idx, bl2] : itertools::enumerate(wdata.gf_struct)) {
39 auto &[bl2_name, bl2_size] = bl2;
40 if (measure_g2w) {
41 g2w[bl1_idx][bl2_idx] = gf<prod<imfreq, imfreq>, tensor_valued<4>>(
42 { mesh_bosonic , mesh_fermionic },
43 make_shape(bl1_size, bl1_size, bl2_size, bl2_size));
44 g2w[bl1_idx][bl2_idx]() = 0;
45 }
46 if (measure_g3w) {
47 g3w[bl1_idx][bl2_idx] = gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4>>(
48 { mesh_bosonic , mesh_fermionic , mesh_fermionic },
49 make_shape(bl1_size, bl1_size, bl2_size, bl2_size));
50 g3w[bl1_idx][bl2_idx]() = 0;
51 }
52 }
53 }
54
55 }
56
57 // -------------------------------------
58
59 void four_point::accumulate(double s) {
60
61 LOG("\n ============ MEASURE FOUR-POINT CORRELATION FUNCTIONS ============ \n");
62
64
88
89 Z += s;
90
91 bool const accumulate_g3w = measure_g3w and (n_cycles_done % g3w_interval == 0);
92 ++n_cycles_done;
93 if (not measure_g2w and not accumulate_g3w) return;
94
95 auto Mw = compute_Mw();
96 auto const &nb_blocks = wdata.gf_struct.size();
97
98 if (measure_g2w) {
99 auto nw = compute_nw();
100 for (auto const &b1 : range(nb_blocks)) {
101 for (auto const &b2 : range(nb_blocks)) {
102 auto const &block_shape = g2w[b1][b2].target_shape();
103 for (auto const &c : range(block_shape[2])) {
104 auto col = wdata.block_to_color(b2, c);
105 for (auto const &a : range(block_shape[0])) {
106 for (auto const &b : range(block_shape[1])) {
107 for (auto const &w : mesh_bosonic) {
108 for (auto const &nu1 : mesh_fermionic) {
109 g2w[b1][b2][w, nu1](a, b, c, c) -= s * Mw[b1][-nu1, nu1 + w](a, b) * nw[col][-w];
110 } // nu1
111 } // w
112 } // b
113 } // a
114 } // c
115 } // b2
116 } // b1
117 } // measure_g2w
118
119 if (accumulate_g3w) {
120 Z_g3w += s;
121 for (auto const &b1 : range(nb_blocks)) {
122 for (auto const &b2 : range(nb_blocks)) {
123 auto const &block_shape = g3w[b1][b2].target_shape();
124 for (auto const &a : range(block_shape[0])) {
125 for (auto const &b : range(block_shape[1])) {
126 for (auto const &c : range(block_shape[2])) {
127 for (auto const &d : range(block_shape[3])) {
128 for (auto const &nu1 : mesh_fermionic) {
129 for (auto const &nu2 : mesh_fermionic) {
130 for (auto const &w : mesh_bosonic) {
131 g3w[b1][b2][w, nu1, nu2](a, b, c, d) += s *
132 Mw[b1][-nu1, nu1 + w](a, b) * Mw[b2][-nu2 - w, nu2.value()](c, d);
133 if (b1 == b2)
134 g3w[b1][b2][w, nu1, nu2](a, b, c, d) -= s *
135 Mw[b1][-nu1, nu2.value()](a, d) * Mw[b2][-nu2 - w, nu1 + w](c, b);
136 } // w
137 } // nu2
138 } // nu1
139 } // d
140 } // c
141 } // b
142 } // a
143 } // b2
144 } // b1
145 } // measure_g3w
146
147 }
148
149 // -------------------------------------
150
151 void four_point::collect_results(mpi::communicator const &c) {
152
153 Z = mpi::all_reduce(Z, c);
154 Z_g3w = mpi::all_reduce(Z_g3w, c);
155
156 if (measure_g2w) g2w = mpi::all_reduce(g2w, c);
157 if (measure_g3w) g3w = mpi::all_reduce(g3w, c);
158 for (auto const &b1 : range(wdata.gf_struct.size())) {
159 for (auto const &b2 : range(wdata.gf_struct.size())) {
160 if (measure_g2w) g2w[b1][b2] = g2w[b1][b2] / (Z * beta);
161 if (measure_g3w) g3w[b1][b2] = g3w[b1][b2] / (Z_g3w * beta);
162 }
163 }
164
165 if (measure_g2w) {
166 g2w_block = make_block2_gf(block_names, block_names, g2w);
167 results.g2w = std::move(g2w_block);
168 }
169 if (measure_g3w) {
170 g3w_block = make_block2_gf(block_names, block_names, g3w);
171 results.g3w = std::move(g3w_block);
172 }
173
174 }
175
176 // -------------------------------------
177
178 block_gf<prod<imfreq, imfreq>> four_point::compute_Mw() {
179
180 // Mw[bl][nu1, nu2](a, b) = \int d \tau_a \, d \tau_b e^{i(\nu_1 \tau_a + \nu_2 \tau_b)} M[bl](\tau_a, \tau_b)
181 // The time tau_a (tau_b) corresponds to a c_dag (c) operator.
182
183 int n_w_aux = n_w_fermionic + n_w_bosonic - 1;
184 auto aux_mesh = triqs::mesh::imfreq(beta, Fermion, n_w_aux);
185 auto Mw = block_gf(prod(aux_mesh, aux_mesh), wdata.gf_struct);
186 Mw() = 0;
187
188 auto w0 = aux_mesh[0].value();
189 auto dw = aux_mesh[1].value() - aux_mesh[0].value();
190
191 for (auto const &[bl, det] : itertools::enumerate(wdata.dets)) {
192 long N = det.size();
193 for (long i : range(N)) {
194 auto [tau_i, a] = det.get_x(i);
195 for (long j : range(N)) {
196 auto [tau_j, b] = det.get_y(j);
197 auto Mij = det.inverse_matrix(j, i); // The row index is i and the column index is j.
198 auto exp_i = std::exp(w0 * double(tau_i));
199 auto exp_i_dw = std::exp(dw * double(tau_i));
200 auto exp_j0 = std::exp(w0 * double(tau_j));
201 auto exp_j_dw = std::exp(dw * double(tau_j));
202 auto exp_j = exp_j0;
203 for (auto const &nu1 : aux_mesh) {
204 exp_j = exp_j0;
205 auto expMij = exp_i * Mij;
206 for (auto const &nu2 : aux_mesh) {
207 Mw[bl][nu1, nu2](a, b) += expMij * exp_j;
208 exp_j = exp_j * exp_j_dw;
209 } // nu2
210 exp_i = exp_i * exp_i_dw;
211 } // nu1
212 } // tau_j
213 } // tau_i
214 } // bl
215 return Mw;
216
217 }
218
219 // -------------------------------------
220
221 std::vector<gf<imfreq, scalar_valued>> four_point::compute_nw() {
222
223 /*
224 $$ n(col)[\omega] = \sum_{\text{segments}}
225 \int_{\tau_\mathrm{start}}^{\tau_\mathrm{end}}e^{i\omega\tau}\mathrm{d}\tau $$
226 */
227
228 std::vector<gf<imfreq, scalar_valued>> nw;
229 nw.resize(wdata.n_color);
230 for (auto const &c : range(wdata.n_color)) {
231 nw[c] = gf<imfreq, scalar_valued>(mesh_bosonic);
232 nw[c]() = 0;
233 nw[c][0] = -beta; // Take care of the convention: nw[omega = 0] = density - 1.
234
235 for (auto const &s: config.seglists[c]) {
236
237 double tau_c = double(s.tau_c);
238 double tau_cdag = double(s.tau_cdag);
239
240 // Zero frequency: Add up the all the segment lengths
241 nw[c][0] += double(s.length());
242
243 // Compute remaining frequencies
244 for (auto const &w : mesh_bosonic) {
245 if (w.n == 0) continue;
246 nw[c][w] += (std::exp(w.value() * tau_c) - std::exp(w.value() * tau_cdag)) / w.value();
247 } // w
248
249 } // s
250 } // c
251 return nw;
252
253 }
254
255} // namespace triqs_ctseg::measures