add four-point vertex measurement

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"

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<gf<prod<imfreq, imfreq, imfreq>, 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<prod<imfreq, imfreq, imfreq>, 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<std::string> block_names;
+    for (auto const &[bl_name, bl_size] : p.gf_struct) block_names.push_back(bl_name);
+    auto bosonic_block_names = std::vector<std::string>{};
+    for (auto const &str1 : block_names)
+      for (auto const &str2 : block_names)
+        bosonic_block_names.push_back(str1 + "|" + str2);
+
+    g3w = make_block_gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4>>(bosonic_block_names, g3w_vec());
+    f3w = make_block_gf<prod<imfreq, imfreq, imfreq>, 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<tau_t, long> 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<array<dcomplex, 4>> 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<array<dcomplex, 4>> 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

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<array<dcomplex, 4>> Mw_vector;
+    vector<array<dcomplex, 4>> nMw_vector;
+    vector<dcomplex> y_exp_ini, y_exp_inc, x_exp_ini, x_exp_inc;
+    vector<int> y_inner_index, x_inner_index;
+
+    block_gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4>> g3w;
+    block_gf<prod<imfreq, imfreq, imfreq>, tensor_valued<4>> f3w;
+
+    double Z = 0;
+
+    four_point(params_t const &params, 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<tau_t, long> const &y);
+
+    vector<array<dcomplex, 4>> 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

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);

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

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