|
| 1 | +// Copyright (c) 2024 Simons Foundation |
| 2 | +// |
| 3 | +// This program is free software: you can redistribute it and/or modify |
| 4 | +// it under the terms of the GNU General Public License as published by |
| 5 | +// the Free Software Foundation, either version 3 of the License, or |
| 6 | +// (at your option) any later version. |
| 7 | +// |
| 8 | +// This program is distributed in the hope that it will be useful, |
| 9 | +// but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 10 | +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 11 | +// GNU General Public License for more details. |
| 12 | +// |
| 13 | +// You may obtain a copy of the License at |
| 14 | +// https://www.gnu.org/licenses/gpl-3.0.txt |
| 15 | +// |
| 16 | +// Authors: Hao Lu |
| 17 | + |
| 18 | +#include "./four_point.hpp" |
| 19 | +#include "../logs.hpp" |
| 20 | + |
| 21 | +namespace triqs_ctseg::measures { |
| 22 | + |
| 23 | + four_point::four_point(params_t const &p, work_data_t const &wdata, configuration_t const &config, results_t &results) |
| 24 | + : wdata{wdata}, config{config}, results{results} { |
| 25 | + |
| 26 | + beta = p.beta; |
| 27 | + measure_g3w = p.measure_g3w; |
| 28 | + measure_f3w = p.measure_f3w; |
| 29 | + |
| 30 | + auto g3w_vec = [&]() { |
| 31 | + std::vector<gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4>>> green_v; |
| 32 | + for (auto const &[bl1_name, bl1_size] : p.gf_struct) |
| 33 | + for (auto const &[bl2_name, bl2_size] : p.gf_struct) |
| 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}}, |
| 38 | + make_shape(bl1_size, bl1_size, bl2_size, bl2_size))); |
| 39 | + return green_v; |
| 40 | + }; |
| 41 | + |
| 42 | + std::vector<std::string> block_names; |
| 43 | + for (auto const &[bl_name, bl_size] : p.gf_struct) block_names.push_back(bl_name); |
| 44 | + auto bosonic_block_names = std::vector<std::string>{}; |
| 45 | + for (auto const &str1 : block_names) |
| 46 | + for (auto const &str2 : block_names) |
| 47 | + bosonic_block_names.push_back(str1 + "|" + str2); |
| 48 | + |
| 49 | + g3w = make_block_gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4>>(bosonic_block_names, g3w_vec()); |
| 50 | + f3w = make_block_gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4>>(bosonic_block_names, g3w_vec()); |
| 51 | + |
| 52 | + for (auto &g : g3w) |
| 53 | + g() = 0; |
| 54 | + for (auto &f : f3w) |
| 55 | + f() = 0; |
| 56 | + |
| 57 | + LOG("\n ====================== COMPUTING M-MATRIX ====================== \n"); |
| 58 | + if (measure_g3w) Mw_vector = compute_Mw(false); |
| 59 | + if (measure_f3w) nMw_vector = compute_Mw(true); |
| 60 | + |
| 61 | + n_w_fermionic = std::get<0>(g3w[0].mesh().components()).last_index() + 1; |
| 62 | + n_w_bosonic = std::get<2>(g3w[0].mesh().components()).last_index() + 1; |
| 63 | + |
| 64 | + w_ini = (2 * (-n_w_fermionic) + 1) * M_PI / beta; |
| 65 | + w_inc = 2 * M_PI / beta; |
| 66 | + |
| 67 | + } |
| 68 | + |
| 69 | + // ------------------------------------- |
| 70 | + |
| 71 | + void four_point::accumulate(double s) { |
| 72 | + |
| 73 | + LOG("\n ============ MEASURE FOUR-POINT CORRELATION FUNCTION ============ \n"); |
| 74 | + |
| 75 | + /// Measure the four-point correlation function |
| 76 | + /** |
| 77 | + *The four-point correlation function is defined as: |
| 78 | + * |
| 79 | + *$$\chi^{\sigma\sigma'}_{abcd}(i\omega, i\omega',i\Omega) = G^{2,\sigma, |
| 80 | + *\sigma'}_{abcd}(i\omega, i\omega',i\Omega) = \langle c_{a\sigma}(i\omega) |
| 81 | + *c^\dagger_{b\sigma}(i\omega+i\Omega) c_{c\sigma'}(i\omega'+i\Omega) |
| 82 | + *c^\dagger_{d\sigma'}(i\omega') \rangle$$ |
| 83 | + * |
| 84 | + * Its improved estimator is the Fourier transform of |
| 85 | + * |
| 86 | + *$$F^{3,\sigma\sigma'}_{abcd}(\tau,\tau',\tau'') = \int \mathrm{d}\bar{\tau} |
| 87 | + *\sum_{e\bar{\sigma}} \mathcal{U}^{\sigma\bar{\sigma}}_{ae}(\bar{\tau}-\tau) |
| 88 | + *\langle c_{a\sigma}(\tau) c^\dagger_{b\sigma}(\tau') c_{c\sigma'}(\tau'') |
| 89 | + *c^\dagger_{d\sigma'}(0) \rangle$$ |
| 90 | + * |
| 91 | + * The vertex corresponding to this correlation function is evaluated separately. |
| 92 | + */ |
| 93 | + |
| 94 | + Z += s; |
| 95 | + |
| 96 | + if (measure_g3w) { |
| 97 | + for (long bl = 0; bl < g3w.size(); bl++) { // bl : 'upup', 'updn', ... |
| 98 | + long b1 = bl / wdata.gf_struct.size(); |
| 99 | + long b2 = bl % wdata.gf_struct.size(); |
| 100 | + for (int a = 0; a < g3w[bl].target_shape()[0]; a++) { |
| 101 | + for (int b = 0; b < g3w[bl].target_shape()[1]; b++) { |
| 102 | + for (int c = 0; c < g3w[bl].target_shape()[2]; c++) { |
| 103 | + for (int d = 0; d < g3w[bl].target_shape()[3]; d++) { |
| 104 | + for (int n1 = -n_w_fermionic; n1 < n_w_fermionic; n1++) { |
| 105 | + for (int n4 = -n_w_fermionic; n4 < n_w_fermionic; n4++) { |
| 106 | + for (int m = 0; m < n_w_bosonic; m++) { |
| 107 | + int n2 = n1 + m; |
| 108 | + int n3 = n4 + m; |
| 109 | + g3w[bl](n1, n4, m)(a, b, c, d) += s * Mw(b1, a, b, n1, n2) * Mw(b2, c, d, n3, n4); |
| 110 | + if (b1 == b2) |
| 111 | + g3w[bl](n1, n4, m)(a, b, c, d) -= s * Mw(b1, a, d, n1, n4) * Mw(b2, c, b, n3, n2); |
| 112 | + } // m |
| 113 | + } // n4 |
| 114 | + } // n1 |
| 115 | + } // d |
| 116 | + } // c |
| 117 | + } // b |
| 118 | + } // a |
| 119 | + } // bl |
| 120 | + } // measure_g3w |
| 121 | + |
| 122 | + if (measure_f3w) { |
| 123 | + for (long bl = 0; bl < f3w.size(); bl++) { // bl : 'upup', 'updn', ... |
| 124 | + long b1 = bl / wdata.gf_struct.size(); |
| 125 | + long b2 = bl % wdata.gf_struct.size(); |
| 126 | + for (int a = 0; a < f3w[bl].target_shape()[0]; a++) { |
| 127 | + for (int b = 0; b < f3w[bl].target_shape()[1]; b++) { |
| 128 | + for (int c = 0; c < f3w[bl].target_shape()[2]; c++) { |
| 129 | + for (int d = 0; d < f3w[bl].target_shape()[3]; d++) { |
| 130 | + for (int n1 = -n_w_fermionic; n1 < n_w_fermionic; n1++) { |
| 131 | + for (int n4 = -n_w_fermionic; n4 < n_w_fermionic; n4++) { |
| 132 | + for (int m = 0; m < n_w_bosonic; m++) { |
| 133 | + int n2 = n1 + m; |
| 134 | + int n3 = n4 + m; |
| 135 | + f3w[bl](n1, n4, m)(a, b, c, d) += s * nMw(b1, a, b, n1, n2) * Mw(b2, c, d, n3, n4); |
| 136 | + if (b1 == b2) |
| 137 | + f3w[bl](n1, n4, m)(a, b, c, d) -= s * nMw(b1, a, d, n1, n4) * Mw(b2, c, b, n3, n2); |
| 138 | + } // m |
| 139 | + } // n4 |
| 140 | + } // n1 |
| 141 | + } // d |
| 142 | + } // c |
| 143 | + } // b |
| 144 | + } // a |
| 145 | + } // bl |
| 146 | + } |
| 147 | + |
| 148 | + } |
| 149 | + |
| 150 | + // ------------------------------------- |
| 151 | + |
| 152 | + void four_point::collect_results(mpi::communicator const &c) { |
| 153 | + |
| 154 | + Z = mpi::all_reduce(Z, c); |
| 155 | + if (measure_g3w) { |
| 156 | + g3w = mpi::all_reduce(g3w, c); |
| 157 | + g3w = g3w / (Z * beta); |
| 158 | + results.g3w = std::move(g3w); |
| 159 | + } |
| 160 | + if (measure_f3w) { |
| 161 | + f3w = mpi::all_reduce(f3w, c); |
| 162 | + f3w = f3w / (Z * beta); |
| 163 | + results.f3w = std::move(f3w); |
| 164 | + } |
| 165 | + |
| 166 | + } |
| 167 | + |
| 168 | + // ------------------------------------- |
| 169 | + |
| 170 | + double four_point::fprefactor(long const &block, std::pair<tau_t, long> const &y) { |
| 171 | + |
| 172 | + int color = wdata.block_to_color(block, y.second); |
| 173 | + double I_tau = 0; |
| 174 | + for (auto const &[c, sl] : itertools::enumerate(config.seglists)) { |
| 175 | + auto ntau = n_tau(y.first, sl); // Density to the right of y.first in sl |
| 176 | + if (c != color) I_tau += wdata.U(c, color) * ntau; |
| 177 | + if (wdata.has_Dt) { |
| 178 | + I_tau -= K_overlap(sl, y.first, false, wdata.Kprime, c, color); |
| 179 | + if (c == color) I_tau -= 2 * real(wdata.Kprime(0)(c, c)); |
| 180 | + } |
| 181 | + if (wdata.has_Jperp) { |
| 182 | + I_tau -= 4 * real(wdata.Kprime_spin(0)(c, color)) * ntau; |
| 183 | + I_tau -= 2 * K_overlap(sl, y.first, false, wdata.Kprime_spin, c, color); |
| 184 | + } |
| 185 | + } |
| 186 | + return I_tau; |
| 187 | + |
| 188 | + } |
| 189 | + |
| 190 | + // ------------------------------------- |
| 191 | + |
| 192 | + vector<array<dcomplex, 4>> four_point::compute_Mw(bool is_nMw) { |
| 193 | + int n_w_aux = 2 * n_w_fermionic + n_w_bosonic > 1 ? 2 * n_w_fermionic + n_w_bosonic - 1 : 0; |
| 194 | + vector<array<dcomplex, 4>> result; |
| 195 | + result.resize(wdata.gf_struct.size()); |
| 196 | + |
| 197 | + for (auto const &[bl, bl_pair] : itertools::enumerate(wdata.gf_struct)) { |
| 198 | + auto const &[bl_name, bl_size] = bl_pair; |
| 199 | + result[bl].resize(make_shape(bl_size, bl_size, n_w_aux, n_w_aux)); |
| 200 | + result[bl]() = 0; |
| 201 | + |
| 202 | + long N = wdata.dets[bl].size(); |
| 203 | + y_exp_ini.resize(N); |
| 204 | + y_exp_inc.resize(N); |
| 205 | + x_exp_ini.resize(N); |
| 206 | + x_exp_inc.resize(N); |
| 207 | + y_inner_index.resize(N); |
| 208 | + x_inner_index.resize(N); |
| 209 | + |
| 210 | + for (long id : range(N)) { |
| 211 | + auto y = wdata.dets[bl].get_y(id); |
| 212 | + auto x = wdata.dets[bl].get_x(id); |
| 213 | + y_exp_ini(id) = std::exp(dcomplex(0, w_ini * double(std::get<0>(y)))); |
| 214 | + y_exp_inc(id) = std::exp(dcomplex(0, w_inc * double(std::get<0>(y)))); |
| 215 | + x_exp_ini(id) = std::exp(dcomplex(0, -w_ini * double(std::get<0>(x)))); |
| 216 | + x_exp_inc(id) = std::exp(dcomplex(0, -w_inc * double(std::get<0>(x)))); |
| 217 | + y_inner_index(id) = std::get<1>(y); |
| 218 | + x_inner_index(id) = std::get<1>(x); |
| 219 | + } |
| 220 | + |
| 221 | + for (long id_y : range(N)) { |
| 222 | + auto y = wdata.dets[bl].get_y(id_y); |
| 223 | + int yj = y_inner_index(id_y); |
| 224 | + double f_fact = is_nMw ? fprefactor(bl, y) : 1.0; |
| 225 | + |
| 226 | + for (long id_x : range(N)) { |
| 227 | + int xi = x_inner_index(id_x); |
| 228 | + dcomplex y_exp = y_exp_ini(id_y); |
| 229 | + dcomplex x_exp = x_exp_ini(id_x); |
| 230 | + auto Minv = wdata.dets[bl].inverse_matrix(id_y, id_x); |
| 231 | + |
| 232 | + for (int n_1 : range(n_w_aux)) { |
| 233 | + for (int n_2 : range(n_w_aux)) { |
| 234 | + auto val = Minv * y_exp * x_exp; |
| 235 | + result[bl](yj, xi, n_1, n_2) += val * f_fact; |
| 236 | + x_exp *= x_exp_inc(id_x); |
| 237 | + } |
| 238 | + x_exp = x_exp_ini(id_x); |
| 239 | + y_exp *= y_exp_inc(id_y); |
| 240 | + } |
| 241 | + } |
| 242 | + } |
| 243 | + } |
| 244 | + |
| 245 | + return result; |
| 246 | + } |
| 247 | + |
| 248 | +} // namespace triqs_ctseg::measures |
0 commit comments