diff --git a/c++/triqs_ctseg/measures.hpp b/c++/triqs_ctseg/measures.hpp index fb44732..013562f 100644 --- a/c++/triqs_ctseg/measures.hpp +++ b/c++/triqs_ctseg/measures.hpp @@ -25,3 +25,4 @@ #include "./measures/average_sign.hpp" #include "./measures/pert_order.hpp" #include "./measures/state_hist.hpp" +#include "./measures/four_point.hpp" diff --git a/c++/triqs_ctseg/measures/four_point.cpp b/c++/triqs_ctseg/measures/four_point.cpp new file mode 100644 index 0000000..fb51bc7 --- /dev/null +++ b/c++/triqs_ctseg/measures/four_point.cpp @@ -0,0 +1,248 @@ +// Copyright (c) 2024 Simons Foundation +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You may obtain a copy of the License at +// https://www.gnu.org/licenses/gpl-3.0.txt +// +// Authors: Hao Lu + +#include "./four_point.hpp" +#include "../logs.hpp" + +namespace triqs_ctseg::measures { + + four_point::four_point(params_t const &p, work_data_t const &wdata, configuration_t const &config, results_t &results) + : wdata{wdata}, config{config}, results{results} { + + beta = p.beta; + measure_g3w = p.measure_g3w; + measure_f3w = p.measure_f3w; + + auto g3w_vec = [&]() { + std::vector, tensor_valued<4>>> green_v; + for (auto const &[bl1_name, bl1_size] : p.gf_struct) + for (auto const &[bl2_name, bl2_size] : p.gf_struct) + green_v.emplace_back(gf, tensor_valued<4>>( + {{beta, Fermion, p.n_iw_chi4_f, imfreq::option::all_frequencies}, + {beta, Fermion, p.n_iw_chi4_f, imfreq::option::all_frequencies}, + {beta, Boson, p.n_iw_chi4_b, imfreq::option::positive_frequencies_only}}, + make_shape(bl1_size, bl1_size, bl2_size, bl2_size))); + return green_v; + }; + + std::vector block_names; + for (auto const &[bl_name, bl_size] : p.gf_struct) block_names.push_back(bl_name); + auto bosonic_block_names = std::vector{}; + for (auto const &str1 : block_names) + for (auto const &str2 : block_names) + bosonic_block_names.push_back(str1 + "|" + str2); + + g3w = make_block_gf, tensor_valued<4>>(bosonic_block_names, g3w_vec()); + f3w = make_block_gf, tensor_valued<4>>(bosonic_block_names, g3w_vec()); + + for (auto &g : g3w) + g() = 0; + for (auto &f : f3w) + f() = 0; + + LOG("\n ====================== COMPUTING M-MATRIX ====================== \n"); + if (measure_g3w) Mw_vector = compute_Mw(false); + if (measure_f3w) nMw_vector = compute_Mw(true); + + n_w_fermionic = std::get<0>(g3w[0].mesh().components()).last_index() + 1; + n_w_bosonic = std::get<2>(g3w[0].mesh().components()).last_index() + 1; + + w_ini = (2 * (-n_w_fermionic) + 1) * M_PI / beta; + w_inc = 2 * M_PI / beta; + + } + + // ------------------------------------- + + void four_point::accumulate(double s) { + + LOG("\n ============ MEASURE FOUR-POINT CORRELATION FUNCTION ============ \n"); + + /// Measure the four-point correlation function + /** + *The four-point correlation function is defined as: + * + *$$\chi^{\sigma\sigma'}_{abcd}(i\omega, i\omega',i\Omega) = G^{2,\sigma, + *\sigma'}_{abcd}(i\omega, i\omega',i\Omega) = \langle c_{a\sigma}(i\omega) + *c^\dagger_{b\sigma}(i\omega+i\Omega) c_{c\sigma'}(i\omega'+i\Omega) + *c^\dagger_{d\sigma'}(i\omega') \rangle$$ + * + * Its improved estimator is the Fourier transform of + * + *$$F^{3,\sigma\sigma'}_{abcd}(\tau,\tau',\tau'') = \int \mathrm{d}\bar{\tau} + *\sum_{e\bar{\sigma}} \mathcal{U}^{\sigma\bar{\sigma}}_{ae}(\bar{\tau}-\tau) + *\langle c_{a\sigma}(\tau) c^\dagger_{b\sigma}(\tau') c_{c\sigma'}(\tau'') + *c^\dagger_{d\sigma'}(0) \rangle$$ + * + * The vertex corresponding to this correlation function is evaluated separately. + */ + + Z += s; + + if (measure_g3w) { + for (long bl = 0; bl < g3w.size(); bl++) { // bl : 'upup', 'updn', ... + long b1 = bl / wdata.gf_struct.size(); + long b2 = bl % wdata.gf_struct.size(); + for (int a = 0; a < g3w[bl].target_shape()[0]; a++) { + for (int b = 0; b < g3w[bl].target_shape()[1]; b++) { + for (int c = 0; c < g3w[bl].target_shape()[2]; c++) { + for (int d = 0; d < g3w[bl].target_shape()[3]; d++) { + for (int n1 = -n_w_fermionic; n1 < n_w_fermionic; n1++) { + for (int n4 = -n_w_fermionic; n4 < n_w_fermionic; n4++) { + for (int m = 0; m < n_w_bosonic; m++) { + int n2 = n1 + m; + int n3 = n4 + m; + g3w[bl](n1, n4, m)(a, b, c, d) += s * Mw(b1, a, b, n1, n2) * Mw(b2, c, d, n3, n4); + if (b1 == b2) + g3w[bl](n1, n4, m)(a, b, c, d) -= s * Mw(b1, a, d, n1, n4) * Mw(b2, c, b, n3, n2); + } // m + } // n4 + } // n1 + } // d + } // c + } // b + } // a + } // bl + } // measure_g3w + + if (measure_f3w) { + for (long bl = 0; bl < f3w.size(); bl++) { // bl : 'upup', 'updn', ... + long b1 = bl / wdata.gf_struct.size(); + long b2 = bl % wdata.gf_struct.size(); + for (int a = 0; a < f3w[bl].target_shape()[0]; a++) { + for (int b = 0; b < f3w[bl].target_shape()[1]; b++) { + for (int c = 0; c < f3w[bl].target_shape()[2]; c++) { + for (int d = 0; d < f3w[bl].target_shape()[3]; d++) { + for (int n1 = -n_w_fermionic; n1 < n_w_fermionic; n1++) { + for (int n4 = -n_w_fermionic; n4 < n_w_fermionic; n4++) { + for (int m = 0; m < n_w_bosonic; m++) { + int n2 = n1 + m; + int n3 = n4 + m; + f3w[bl](n1, n4, m)(a, b, c, d) += s * nMw(b1, a, b, n1, n2) * Mw(b2, c, d, n3, n4); + if (b1 == b2) + f3w[bl](n1, n4, m)(a, b, c, d) -= s * nMw(b1, a, d, n1, n4) * Mw(b2, c, b, n3, n2); + } // m + } // n4 + } // n1 + } // d + } // c + } // b + } // a + } // bl + } + + } + + // ------------------------------------- + + void four_point::collect_results(mpi::communicator const &c) { + + Z = mpi::all_reduce(Z, c); + if (measure_g3w) { + g3w = mpi::all_reduce(g3w, c); + g3w = g3w / (Z * beta); + results.g3w = std::move(g3w); + } + if (measure_f3w) { + f3w = mpi::all_reduce(f3w, c); + f3w = f3w / (Z * beta); + results.f3w = std::move(f3w); + } + + } + + // ------------------------------------- + + double four_point::fprefactor(long const &block, std::pair const &y) { + + int color = wdata.block_to_color(block, y.second); + double I_tau = 0; + for (auto const &[c, sl] : itertools::enumerate(config.seglists)) { + auto ntau = n_tau(y.first, sl); // Density to the right of y.first in sl + if (c != color) I_tau += wdata.U(c, color) * ntau; + if (wdata.has_Dt) { + I_tau -= K_overlap(sl, y.first, false, wdata.Kprime, c, color); + if (c == color) I_tau -= 2 * real(wdata.Kprime(0)(c, c)); + } + if (wdata.has_Jperp) { + I_tau -= 4 * real(wdata.Kprime_spin(0)(c, color)) * ntau; + I_tau -= 2 * K_overlap(sl, y.first, false, wdata.Kprime_spin, c, color); + } + } + return I_tau; + + } + + // ------------------------------------- + + vector> four_point::compute_Mw(bool is_nMw) { + int n_w_aux = 2 * n_w_fermionic + n_w_bosonic > 1 ? 2 * n_w_fermionic + n_w_bosonic - 1 : 0; + vector> result; + result.resize(wdata.gf_struct.size()); + + for (auto const &[bl, bl_pair] : itertools::enumerate(wdata.gf_struct)) { + auto const &[bl_name, bl_size] = bl_pair; + result[bl].resize(make_shape(bl_size, bl_size, n_w_aux, n_w_aux)); + result[bl]() = 0; + + long N = wdata.dets[bl].size(); + y_exp_ini.resize(N); + y_exp_inc.resize(N); + x_exp_ini.resize(N); + x_exp_inc.resize(N); + y_inner_index.resize(N); + x_inner_index.resize(N); + + for (long id : range(N)) { + auto y = wdata.dets[bl].get_y(id); + auto x = wdata.dets[bl].get_x(id); + y_exp_ini(id) = std::exp(dcomplex(0, w_ini * double(std::get<0>(y)))); + y_exp_inc(id) = std::exp(dcomplex(0, w_inc * double(std::get<0>(y)))); + x_exp_ini(id) = std::exp(dcomplex(0, -w_ini * double(std::get<0>(x)))); + x_exp_inc(id) = std::exp(dcomplex(0, -w_inc * double(std::get<0>(x)))); + y_inner_index(id) = std::get<1>(y); + x_inner_index(id) = std::get<1>(x); + } + + for (long id_y : range(N)) { + auto y = wdata.dets[bl].get_y(id_y); + int yj = y_inner_index(id_y); + double f_fact = is_nMw ? fprefactor(bl, y) : 1.0; + + for (long id_x : range(N)) { + int xi = x_inner_index(id_x); + dcomplex y_exp = y_exp_ini(id_y); + dcomplex x_exp = x_exp_ini(id_x); + auto Minv = wdata.dets[bl].inverse_matrix(id_y, id_x); + + for (int n_1 : range(n_w_aux)) { + for (int n_2 : range(n_w_aux)) { + auto val = Minv * y_exp * x_exp; + result[bl](yj, xi, n_1, n_2) += val * f_fact; + x_exp *= x_exp_inc(id_x); + } + x_exp = x_exp_ini(id_x); + y_exp *= y_exp_inc(id_y); + } + } + } + } + + return result; + } + +} // namespace triqs_ctseg::measures diff --git a/c++/triqs_ctseg/measures/four_point.hpp b/c++/triqs_ctseg/measures/four_point.hpp new file mode 100644 index 0000000..d40abb1 --- /dev/null +++ b/c++/triqs_ctseg/measures/four_point.hpp @@ -0,0 +1,64 @@ +// Copyright (c) 2024 Simons Foundation +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You may obtain a copy of the License at +// https://www.gnu.org/licenses/gpl-3.0.txt +// +// Authors: Hao Lu + +#pragma once +#include "../configuration.hpp" +#include "../work_data.hpp" +#include "../results.hpp" + +namespace triqs_ctseg::measures { + + struct four_point { + + work_data_t const &wdata; + configuration_t const &config; + results_t &results; + double beta; + double w_ini, w_inc; + bool measure_g3w; + bool measure_f3w; + int n_w_fermionic; + int n_w_bosonic; + vector> Mw_vector; + vector> nMw_vector; + vector y_exp_ini, y_exp_inc, x_exp_ini, x_exp_inc; + vector y_inner_index, x_inner_index; + + block_gf, tensor_valued<4>> g3w; + block_gf, tensor_valued<4>> f3w; + + double Z = 0; + + four_point(params_t const ¶ms, work_data_t const &wdata, configuration_t const &config, results_t &results); + + void accumulate(double s); + void collect_results(mpi::communicator const &c); + double fprefactor(long const &block, std::pair const &y); + + vector> compute_Mw(bool is_nMw); + + dcomplex Mw(long const &block, int const &i, int const &j, int const &n1, int const &n2) { + return Mw_vector[block](i, j, n1 + n_w_fermionic, n2 + n_w_fermionic); + } + + dcomplex nMw(long const &block, int const &i, int const &j, int const &n1, int const &n2) { + return nMw_vector[block](i, j, n1 + n_w_fermionic, n2 + n_w_fermionic); + } + + }; + +} // namespace triqs_ctseg::measures \ No newline at end of file diff --git a/c++/triqs_ctseg/params.cpp b/c++/triqs_ctseg/params.cpp index b912109..67d5f11 100644 --- a/c++/triqs_ctseg/params.cpp +++ b/c++/triqs_ctseg/params.cpp @@ -51,6 +51,8 @@ namespace triqs_ctseg { h5_write(grp, "h_loc0", c.h_loc0); h5_write(grp, "n_tau_G", c.n_tau_G); h5_write(grp, "n_tau_chi2", c.n_tau_chi2); + h5_write(grp, "n_iw_chi4_b", c.n_iw_chi4_b); + h5_write(grp, "n_iw_chi4_f", c.n_iw_chi4_f); h5_write(grp, "n_cycles", c.n_cycles); h5_write(grp, "length_cycle", c.length_cycle); h5_write(grp, "n_warmup_cycles", c.n_warmup_cycles); @@ -77,6 +79,8 @@ namespace triqs_ctseg { h5_write(grp, "measure_nn_tau", c.measure_nn_tau); h5_write(grp, "measure_Sperp_tau", c.measure_Sperp_tau); h5_write(grp, "measure_state_hist", c.measure_state_hist); + h5_write(grp, "measure_g3w", c.measure_g3w); + h5_write(grp, "measure_f3w", c.measure_f3w); h5_write(grp, "det_init_size", c.det_init_size); h5_write(grp, "det_n_operations_before_check", c.det_n_operations_before_check); h5_write(grp, "det_precision_warning", c.det_precision_warning); @@ -121,6 +125,8 @@ namespace triqs_ctseg { h5_read(grp, "measure_nn_tau", c.measure_nn_tau); h5_read(grp, "measure_Sperp_tau", c.measure_Sperp_tau); h5_read(grp, "measure_state_hist", c.measure_state_hist); + h5_read(grp, "measure_g3w", c.measure_g3w); + h5_read(grp, "measure_f3w", c.measure_f3w); h5_read(grp, "det_init_size", c.det_init_size); h5_read(grp, "det_n_operations_before_check", c.det_n_operations_before_check); h5_read(grp, "det_precision_warning", c.det_precision_warning); diff --git a/c++/triqs_ctseg/params.hpp b/c++/triqs_ctseg/params.hpp index 59b89cf..0566f41 100644 --- a/c++/triqs_ctseg/params.hpp +++ b/c++/triqs_ctseg/params.hpp @@ -57,6 +57,12 @@ namespace triqs_ctseg { /// Number of points on which to measure 2-point functions (defaults to n_tau_bosonic) int n_tau_chi2 = 0; + /// Number of bosonic M-frequency points on which to measure 4-point functions + int n_iw_chi4_b = 1001; + + /// Number of fermionic M-frequency points on which to measure 4-point functions + int n_iw_chi4_f = 1000; + /// Number of QMC cycles int n_cycles; @@ -139,6 +145,12 @@ namespace triqs_ctseg { /// Whether to measure state histograms (see measures/state_hist) bool measure_state_hist = false; + /// Whether to measure four-point correlation function (see measures/four_point) + bool measure_g3w = false; + + /// Whether to measure four-point correlation function improved estimator (see measures/four_point) + bool measure_f3w = false; + // -------- Misc parameters -------------- /// The maximum size of the determinant matrix before a resize diff --git a/c++/triqs_ctseg/results.cpp b/c++/triqs_ctseg/results.cpp index 6e9fcba..b53a791 100644 --- a/c++/triqs_ctseg/results.cpp +++ b/c++/triqs_ctseg/results.cpp @@ -35,6 +35,8 @@ namespace triqs_ctseg { h5_write(grp, "pert_order_Jperp", c.pert_order_Jperp); h5_write(grp, "average_order_Jperp", c.average_order_Jperp); h5_write(grp, "state_hist", c.state_hist); + h5_write(grp, "g3w", c.g3w); + h5_write(grp, "f3w", c.f3w); } //------------------------------------ @@ -55,6 +57,8 @@ namespace triqs_ctseg { h5_read(grp, "pert_order_Jperp", c.pert_order_Jperp); h5_read(grp, "average_order_Jperp", c.average_order_Jperp); h5_read(grp, "state_hist", c.state_hist); + h5_read(grp, "g3w", c.g3w); + h5_read(grp, "f3w", c.f3w); } } // namespace triqs_ctseg diff --git a/c++/triqs_ctseg/results.hpp b/c++/triqs_ctseg/results.hpp index 0bfaa47..edd0d23 100644 --- a/c++/triqs_ctseg/results.hpp +++ b/c++/triqs_ctseg/results.hpp @@ -60,6 +60,12 @@ namespace triqs_ctseg { /// State histogram std::optional> state_hist; + /// Four-point correlation function + std::optional, tensor_valued<4>>> g3w; + + /// Four-point correlation function improved estimator + std::optional, tensor_valued<4>>> f3w; + /// Average sign double average_sign; }; diff --git a/c++/triqs_ctseg/solver_core.cpp b/c++/triqs_ctseg/solver_core.cpp index 52d889d..f0ed2de 100644 --- a/c++/triqs_ctseg/solver_core.cpp +++ b/c++/triqs_ctseg/solver_core.cpp @@ -125,6 +125,8 @@ namespace triqs_ctseg { } } if (p.measure_state_hist) CTQMC.add_measure(measures::state_hist{p, wdata, config, results}, "State histograms"); + if (p.measure_g3w || p.measure_f3w) + CTQMC.add_measure(measures::four_point{p, wdata, config, results}, "Four-point correlation function"); // Run and collect results CTQMC.warmup_and_accumulate(p.n_warmup_cycles, p.n_cycles, p.length_cycle, diff --git a/python/triqs_ctseg/solver_core_desc.py b/python/triqs_ctseg/solver_core_desc.py index 4c75d23..38028f0 100644 --- a/python/triqs_ctseg/solver_core_desc.py +++ b/python/triqs_ctseg/solver_core_desc.py @@ -91,10 +91,15 @@ read_only= True, doc = r"""State histogram""") -c.add_member(c_name = "average_sign", - c_type = "double", +c.add_member(c_name = "g3w", + c_type = "std::optional, tensor_valued<4>>>", + read_only= True, + doc = r"""Four-point correlation function""") + +c.add_member(c_name = "f3w", + c_type = "std::optional, tensor_valued<4>>>", read_only= True, - doc = r"""Average sign""") + doc = r"""Four-point correlation function improved estimator""") module.add_class(c) @@ -394,6 +399,16 @@ initializer = """ false """, doc = r"""Whether to measure state histograms (see measures/state_hist)""") +c.add_member(c_name = "measure_g3w", + c_type = "bool", + initializer = """ false """, + doc = r"""Whether to measure four-point correlation function (see measures/four_point)""") + +c.add_member(c_name = "measure_f3w", + c_type = "bool", + initializer = """ false """, + doc = r"""Whether to measure four-point correlation function improved estimator (see measures/four_point)""") + c.add_member(c_name = "det_init_size", c_type = "int", initializer = """ 100 """,