[measure] Add direct measure for G_iw

.clang-tidy

@@ -1,2 +1,6 @@
+<<<<<<< HEAD
 Checks: '-*,modernize-*,cppcoreguidelines-*,bugprone-*,-modernize-use-trailing-return-type,-cppcoreguidelines-special-member-functions,-cppcoreguidelines-macro-usage,-cppcoreguidelines-no-malloc,-cppcoreguidelines-pro-bounds-pointer-arithmetic,-cppcoreguidelines-pro-bounds-constant-array-index,-cppcoreguidelines-avoid-magic-numbers, -cppcoreguidelines-non-private-member-variables-in-classes, -bugprone-easily-swappable-parameters'
+=======
+Checks: '-*,modernize-*,cppcoreguidelines-*,-modernize-use-trailing-return-type,-cppcoreguidelines-avoid-magic-numbers,-cppcoreguidelines-virtual-class-destructor'
+>>>>>>> ccac36ea ([measure] Add direct measure for G_iw)
 HeaderFilterRegex: 'triqs_cthyb'

.gitignore

@@ -1,2 +1,8 @@
 compile_commands.json
 doc/cpp2rst_generated
+.*.sw?
+code.vim
+.cache
+.vscode
+.clangd
+

c++/triqs_cthyb/container_set.cpp

@@ -31,6 +31,7 @@ namespace triqs_cthyb {
     h5_write(grp, "G_tau", c.G_tau);
     h5_write(grp, "G_tau_accum", c.G_tau_accum);
     h5_write(grp, "asymmetry_G_tau", c.asymmetry_G_tau);
+    h5_write(grp, "G_iw_direct", c.G_iw_direct);
     h5_write(grp, "G_l", c.G_l);
     h5_write(grp, "O_tau", c.O_tau);
     h5_write(grp, "perturbation_order", c.perturbation_order);
@@ -56,6 +57,7 @@ namespace triqs_cthyb {
     h5_read(grp, "G_tau_accum", c.G_tau_accum);
     h5_read(grp, "asymmetry_G_tau", c.asymmetry_G_tau);
     h5_read(grp, "G_l", c.G_l);
+    h5_read(grp, "G_iw_direct", c.G_iw_direct);
     h5::try_read(grp, "O_tau", c.O_tau);
     h5::try_read(grp, "perturbation_order", c.perturbation_order);
     h5::try_read(grp, "perturbation_order_total", c.perturbation_order_total);

c++/triqs_cthyb/container_set.hpp

@@ -29,7 +29,7 @@ namespace triqs_cthyb {
 
   // Containers for measurements
   struct container_set_t {
-    
+
     // -- Single particle Green's functions
 
     /// Single-particle Green's function :math:`G(\tau)` in imaginary time.
@@ -41,6 +41,9 @@ namespace triqs_cthyb {
     /// Violation of the fundamental Green function property G(tau)[i,j] = G(tau)*[j,i] after the measurement
     std::optional<G_tau_G_target_t> asymmetry_G_tau;
 
+    /// Single-particle Green's function :math:`G(\tau)` in imaginary time.
+    std::optional<G_iw_t> G_iw_direct;
+
     /// Single-particle Green's function :math:`G_l` in Legendre polynomial representation.
     std::optional<G_l_t> G_l;
 

c++/triqs_cthyb/measures/G_iw.cpp

@@ -0,0 +1,68 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2014, H. U.R. Strand, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet
+ *
+ * TRIQS 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.
+ *
+ * TRIQS 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 should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "./G_iw.hpp"
+#include <triqs/mesh/imfreq.hpp>
+
+namespace triqs_cthyb {
+
+  using namespace triqs::gfs;
+  using namespace triqs::mesh;
+
+  measure_G_iw ::measure_G_iw(qmc_data const &data, int n_iw, gf_struct_t const &gf_struct, container_set_t &results) : data(data), average_sign(0) {
+    results.G_iw_direct = block_gf<imfreq, G_target_t>{{data.config.beta(), Fermion, n_iw}, gf_struct};
+    G_iw.rebind(*results.G_iw_direct);
+    G_iw() = 0.0;
+  }
+
+  void measure_G_iw::accumulate(mc_weight_t s) {
+    s *= data.atomic_reweighting;
+    average_sign += s;
+
+    for (auto block_idx : range(G_iw.size())) {
+      foreach (data.dets[block_idx], [this, s, block_idx](op_t const &x, op_t const &y, det_scalar_t M) {
+        // beta-periodicity is implicit in the argument, just fix the sign properly
+        auto val      = (y.first >= x.first ? s : -s) * M;
+        double dtau   = double(y.first - x.first);
+        auto const &m = G_iw[0].mesh();
+        long n_iw     = m.last_index();
+        auto ex_dt    = std::exp(M_PI / m.domain().beta * dtau * 1i);
+        auto ex_2dt   = ex_dt * ex_dt;
+        auto ex       = ex_dt;
+        for (long n = 0; n < n_iw; ++n, ex *= ex_2dt) { G_iw[block_idx][n](y.second, x.second) += val * ex; }
+        //this->G_tau[block_idx][closest_mesh_pt(dtau)](y.second, x.second) += val;
+      })
+        ;
+    }
+  }
+
+  void measure_G_iw::collect_results(mpi::communicator const &c) {
+
+    G_iw         = mpi::all_reduce(G_iw, c);
+    average_sign = mpi::all_reduce(average_sign, c);
+    G_iw /= -real(average_sign);
+    // for (auto &gbl : G_iw) {
+    //   //double beta = gbl.mesh().domain().beta;
+    //   gbl /= -real(average_sign); //  * beta * gbl.mesh().delta();
+    // }
+  }
+
+} // namespace triqs_cthyb

c++/triqs_cthyb/measures/G_iw.hpp

@@ -0,0 +1,50 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2014, H. U.R. Strand, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet
+ *
+ * TRIQS 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.
+ *
+ * TRIQS 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 should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+
+#include <triqs/gfs.hpp>
+#include <triqs/mesh.hpp>
+
+#include "../qmc_data.hpp"
+#include "../container_set.hpp"
+
+namespace triqs_cthyb {
+
+  using namespace triqs::gfs;
+  using namespace triqs::mesh;
+
+  // Measure imaginary time Green's function (all blocks)
+  class measure_G_iw {
+
+    public:
+    measure_G_iw(qmc_data const &data, int n_iw, gf_struct_t const &gf_struct, container_set_t &results);
+    void accumulate(mc_weight_t s);
+    void collect_results(mpi::communicator const &c);
+
+    private:
+    qmc_data const &data; //NOLINT
+    mc_weight_t average_sign;
+    G_iw_t::view_type G_iw;
+  };
+
+  // FIXME : assert concept
+
+} // namespace triqs_cthyb

c++/triqs_cthyb/measures/G_tau.cpp

@@ -61,12 +61,12 @@ namespace triqs_cthyb {
       G_tau_block /= -real(average_sign) * beta * G_tau_block.mesh().delta();
 
       // Multiply first and last bins by 2 to account for full bins
-      int last = G_tau_block.mesh().size() - 1;
+      long last = G_tau_block.mesh().size() - 1;
       G_tau_block[0] *= 2;
       G_tau_block[last] *= 2;
 
       // Enforce discontinuity in Green function
-      G_tau_block[0] = 0.5 * matrix_t(G_tau_block[0] - 1 - G_tau_block[last]);
+      G_tau_block[0]    = 0.5 * matrix_t(G_tau_block[0] - 1 - G_tau_block[last]);
       G_tau_block[last] = -1 - G_tau_block[0];
     }
 

c++/triqs_cthyb/measures/G_tau.hpp

@@ -40,7 +40,7 @@ namespace triqs_cthyb {
     void collect_results(mpi::communicator const &c);
 
     private:
-    qmc_data const &data;
+    qmc_data const &data; //NOLINT
     mc_weight_t average_sign;
     G_tau_G_target_t::view_type G_tau;
     G_tau_G_target_t::view_type asymmetry_G_tau;

c++/triqs_cthyb/parameters.hpp

@@ -122,6 +122,12 @@ namespace triqs_cthyb {
     /// Measure G(tau)? :math:`G_{ij}(\tau)=G_{ji}^*(\tau)` is enforced for the resulting G(tau)
     bool measure_G_tau = true;
 
+    /// Measure G(iw) directly
+    bool measure_G_iw = false;
+
+    /// Number of Matsubara frequencies in the direct measurement of G(iw)
+    int measure_G_n_iw = 20;
+
     /// Measure G_l (Legendre)? Note, no hermiticity in G_l is ensured
     bool measure_G_l = false;
 
@@ -248,6 +254,6 @@ namespace triqs_cthyb {
   /// A struct combining both constr_params_t and solve_params_t
   struct params_t : constr_parameters_t, solve_parameters_t {
     params_t(constr_parameters_t constr_parameters_, solve_parameters_t solve_parameters_)
-      : constr_parameters_t(constr_parameters_), solve_parameters_t(solve_parameters_) {}
+       : constr_parameters_t(constr_parameters_), solve_parameters_t(solve_parameters_) {}
   };
-}
+} // namespace triqs_cthyb