Program Listing for File average.cpp

↰ Return to documentation for file (src_qcm/average.cpp)

#include <fstream>
#include "lattice_model_instance.hpp"
#include "integrate.hpp"
#include "parser.hpp"
#include "QCM.hpp"

double tr_sigma_inf(0.0);

vector<shared_ptr<lattice_operator>> ops = {};

extern vector<double> grid_freqs; // optional imaginary frequency grid for discrete integrals
extern vector<double> grid_weights; // weights associated with grid_freqs

//==============================================================================
vector<pair<string,double>> lattice_model_instance::averages(const vector<string> &_ops)
{
  if(average_solved) return ave;
  if(_ops.size() == 0) ops = model->one_body_ops;
  else{
    ops.clear();
    ops.reserve(_ops.size());
    for(int i=0; i<_ops.size(); i++){
      if(model->term[_ops[i]]->is_interaction == false) ops.push_back(model->term[_ops[i]]);
    }
  }
  if(!gf_solved) Green_function_solve();
  double accur_OP = global_double("accur_OP");
  bool periodized_averages = global_bool("periodized_averages");


    vector<double> Iv(ops.size());
  vector<double> Iw(ops.size());

  if(model->hybrid != nullptr){ // frequency-momentum sum
      vector<double> I(ops.size());
    for(int iw=0; iw < model->hybrid->nw; iw++){
      to_zero(Iw);
      for(int ik=0; ik < model->hybrid->nk; ik++){
        average_integrand(iw, ik, I);
        Iw += I;
      }
      Iw *= model->hybrid->weight[iw]/model->hybrid->nk;
      Iv += Iw;
    }
  }
    else{ // frequency-momentum integral
    // lambda function
    if(periodized_averages){
      auto F = [this] (Complex w, vector3D<double> &k, const int *nv, double *I) mutable {average_integrand_per(w, k, nv, I);};
      QCM::wk_integral(model->spatial_dimension, F, Iv, accur_OP, global_bool("verb_integrals"));
    }
    else if(grid_freqs.size() > 0){
      // Frequency-momentum grid. The cluster Green function depends only on w,
      // so build it once per frequency (sequential), then parallelize only over k
      // via k_integral_grid. The lambda captures G_up (and G_dn) by value: the
      // captured copies live in the lambda object, are read-only inside the parallel
      // section, and contain no shared mutable state.
      int nkx, nky, nkz;
      QCM::get_wavevector_grid(nkx, nky, nkz);
      cout << "computing integrals from a grid of " << grid_freqs.size() << " frequencies and "
           << nkx << "x" << nky << "x" << nkz << " k-points (dim=" << model->spatial_dimension << ")" << endl;
      bool has_dn = (model->mixing == HS_mixing::up_down);
      for(int iw = 0; iw < (int)grid_freqs.size(); iw++){
        Complex wc(0, grid_freqs[iw]);
        Green_function G_up = cluster_Green_function(wc, false, false);
        Green_function G_dn;
        if(has_dn) G_dn = cluster_Green_function(wc, false, true);

        auto F_k = [this, G_up, G_dn, has_dn] (vector3D<double> &k, const int *nv, double *I) mutable {
          average_integrand_k(G_up, has_dn ? &G_dn : nullptr, k, nv, I);
        };
        to_zero(Iw);
        QCM::k_integral_grid(model->spatial_dimension, nkx, nky, nkz, F_k, Iw);
        Iv += Iw*grid_weights[iw];
      }
    }
    else{
      auto F = [this] (Complex w, vector3D<double> &k, const int *nv, double *I) mutable {average_integrand(w, k, nv, I);};
      QCM::wk_integral(model->spatial_dimension, F, Iv, accur_OP, global_bool("verb_integrals"));
    }
  }


  size_t i = 0;
  ave.resize(ops.size());
    for(auto& op : ops){
    if(periodized_averages) Iv[i] *= model->Lc;
    if(model->mixing == HS_mixing::full) {
      Iv[i] += op->nambu_correction_full;
      Iv[i] *= 0.5;
    }
    else if((model->mixing&HS_mixing::anomalous)  == HS_mixing::anomalous) Iv[i] += op->nambu_correction;
    ave[i] = {op->name, Iv[i]*op->norm};
        i++;
    }
  average_solved = true;

  // computing the kinetic energy
  E_kin = 0.0;
  i = 0;
  for(auto& op : ops){
    if(op->name != "mu") E_kin += Iv[i]*params.at(op->name);
    i++;
  }
  E_kin /= model->Lc;

    return ave;
}




//==============================================================================
void lattice_model_instance::average_integrand(Complex w, vector3D<double> &k, const int *nv, double *I)
{
  check_signals();
  Green_function G_up = cluster_Green_function(w, false, false);
  if(model->mixing == HS_mixing::up_down){
    Green_function G_dn = cluster_Green_function(w, false, true);
    average_integrand_k(G_up, &G_dn, k, nv, I);
  }
  else average_integrand_k(G_up, nullptr, k, nv, I);
}


//==============================================================================
void lattice_model_instance::average_integrand_k(Green_function &G_up, Green_function *G_down, vector3D<double> &k, const int *nv, double *I)
{
  Complex w = G_up.w;
  const double w_offset = 2.0;
  Complex G_pole = 1.0/(w - w_offset);

  Green_function_k K(G_up,k);
  set_Gcpt(K);
  K.Gcpt.add(-G_pole); // regulates the Green function at high frequency

  size_t i = 0;
  for(auto& op : ops){
    Complex z(0.0);
    for(auto &x : op->GF_elem) z += K.Gcpt(x.c,x.r)*x.v;
    for(auto &x : op->IGF_elem) z += K.Gcpt(x.c,x.r)*x.v*K.phase[x.n];
    I[i++] = real<double>(z);
  }

  if(G_down){
    Green_function_k K_dn(*G_down,k);
    set_Gcpt(K_dn);
    K_dn.Gcpt.add(-G_pole);
    i = 0;
    for(auto& op : ops){
      Complex z(0.0);
      for(auto &x : op->GF_elem_down) z += K_dn.Gcpt(x.c,x.r)*x.v;
      for(auto &x : op->IGF_elem_down) z += K_dn.Gcpt(x.c,x.r)*x.v*K_dn.phase[x.n];
      I[i++] += real<double>(z);
    }
  }
  if(model->mixing == HS_mixing::normal) for(size_t i=0; i<ops.size(); i++) I[i] *= 2;
}



//==============================================================================
void lattice_model_instance::average_integrand(int iw, int ik, vector<double> &I)
{

  Complex w(0,model->hybrid->w[iw]);
  const double w_offset = 2.0;
  Complex G_pole = 1.0/(w - w_offset);

  Green_function G = cluster_Green_function(w, false, false);
  G.iw = iw;
    Green_function_k K(G, model->hybrid->k[ik], ik);
    set_Gcpt(K);
    K.Gcpt.add(-G_pole); // regulates the Green function at high frequency (subtracts G_pole times the identity matrix)

    size_t i = 0;
    for(auto& op : ops){
        Complex z(0.0);
        for(auto &x : op->GF_elem) z += K.Gcpt(x.c,x.r)*x.v;
        for(auto &x : op->IGF_elem) z += K.Gcpt(x.c,x.r)*x.v*K.phase[x.n];

        I[i++] = real<double>(z);
    }
    i = 0;
    if(model->mixing == HS_mixing::up_down){
    Green_function G = cluster_Green_function(w, false, true);
    Green_function_k K(G, model->hybrid->k[ik], ik);
    set_Gcpt(K);
    K.Gcpt.add(-G_pole); // regulates the Green function at high frequency (subtracts G_pole times the identity matrix)
        for(auto& op : ops){
            Complex z(0.0);
            for(auto &x : op->GF_elem_down) z += K.Gcpt(x.c,x.r)*x.v;
            for(auto &x : op->IGF_elem_down) z += K.Gcpt(x.c,x.r)*x.v*K.phase[x.n];
            I[i++] += real<double>(z);
        }
    }
  if(model->mixing == HS_mixing::normal) for(size_t i=0; i<ops.size(); i++) I[i] *= 2;
}

//==============================================================================
void lattice_model_instance::average_integrand_per(Complex w, vector3D<double> &k, const int *nv, double *I)
{

  const double w_offset = 2.0;
  Complex G_pole = 1.0/(w - w_offset);

  Green_function G = cluster_Green_function(w, false, false);
  vector3D<double> k_red = model->superdual.to(model->physdual.from(k));

  Green_function_k K(G,k_red);
  periodized_Green_function(K);
  K.g.add(-G_pole); // regulates the Green function at high frequency (subtracts G_pole times the identity matrix)

  size_t i = 0;
  for(auto& op : ops){
    matrix<Complex> S(G.G.r);
    // build the operator matrix
    for(auto &x : op->GF_elem) S(x.r, x.c) += x.v;
    for(auto &x : op->IGF_elem) S(x.r, x.c) += x.v*K.phase[x.n];
    matrix<Complex> S_per = model->periodize(k_red, S);
    I[i++] = realpart(K.g.trace_product(S_per));
  }
  i = 0;
  if(model->mixing == HS_mixing::up_down){
    Green_function G = cluster_Green_function(w,false,true);
    Green_function_k K(G,k_red);
    periodized_Green_function(K);
    K.g.add(-G_pole); // regulates the Green function at high frequency (subtracts G_pole times the identity matrix)
    for(auto& op : ops){
      matrix<Complex> S(G.G.r);
      // build the operator matrix
      for(auto &x : op->GF_elem_down) S(x.r, x.c) += x.v;
      for(auto &x : op->IGF_elem_down) S(x.r, x.c) += x.v*K.phase[x.n];
      matrix<Complex> S_per = model->periodize(k_red, S);
      I[i++] += realpart(K.g.trace_product(S_per));
    }
  }
  if(model->mixing == HS_mixing::normal) for(size_t i=0; i<ops.size(); i++) I[i] *= 2;
}





//==============================================================================
vector<double> lattice_model_instance::dos(const complex<double> w, bool use_grid)
{
  double accur_OP = global_double("accur_OP");
  if(!gf_solved) Green_function_solve();

  size_t D_dim = model->n_band;
  size_t d = D_dim;
  if(model->mixing&3) d *= 2;
  vector<double> D(2*D_dim);

  Green_function G = cluster_Green_function(w, false, false);

  auto F = [this, G] (vector3D<double> &k, const int *nv, double *I) mutable {
    for(size_t i = 0 ; i<*nv; i++) I[i] = 0.0;
    Green_function_k K(G,k);
    set_Gcpt(K);
    for(size_t i=0; i<model->dim_GF; i++) I[model->reduced_Green_index[i]] += K.Gcpt(i,i).imag();
  };

  vector<double> Iv(model->dim_reduced_GF,0.0);
  if(use_grid){
    int nkx, nky, nkz;
    QCM::get_wavevector_grid(nkx, nky, nkz);
    QCM::k_integral_grid(model->spatial_dimension, nkx, nky, nkz, F, Iv);
  }
  else QCM::k_integral(model->spatial_dimension, F, Iv, accur_OP, global_bool("verb_integrals"));
  for(size_t i=0; i<d; i++) D[i] = -M_1_PI*Iv[i]/model->Lc;

  if(model->mixing == HS_mixing::up_down){
    Green_function G_down = cluster_Green_function(w, false, true);

    auto F_down = [this, G_down] (vector3D<double> &k, const int *nv, double *I) mutable {
      for(size_t i = 0 ; i<*nv; i++) I[i] = 0.0;
      Green_function_k K(G_down,k);
      set_Gcpt(K);
      for(size_t i=0; i<model->dim_GF; i++) I[model->reduced_Green_index[i]] += K.Gcpt(i,i).imag();
    };

    to_zero(Iv);
    if(use_grid){
      int nkx, nky, nkz;
      QCM::get_wavevector_grid(nkx, nky, nkz);
      QCM::k_integral_grid(model->spatial_dimension, nkx, nky, nkz, F_down, Iv);
    }
    else QCM::k_integral(model->spatial_dimension, F_down, Iv, accur_OP, global_bool("verb_integrals"));
    for(size_t i=0; i<D_dim; i++) D[i+D_dim] = -M_1_PI*Iv[i]/model->Lc;
  }
  else if(model->mixing == HS_mixing::normal){
    for(size_t i=0; i<D_dim; i++) D[i+D_dim] = D[i];
  }
  return D;
}


//==============================================================================
double lattice_model_instance::spectral_average(const string& name, const complex<double> w)
{
  double accur_OP = global_double("accur_OP");
  if(model->term.find(name) == model->term.end()) qcm_throw("the operator "+name+" does not exist");
  size_t n_clus = model->clusters.size();

  lattice_operator op = *model->term.at(name);

  if(!gf_solved) Green_function_solve();

  Green_function G = cluster_Green_function(w, false, false);


  auto F = [this, op, G] (vector3D<double> &k, const int *nv, double *I) mutable {
    for(size_t i = 0 ; i<*nv; i++) I[i] = 0.0;
    Green_function_k K(G,k);
    set_Gcpt(K);
    Complex z(0.0);
    for(auto &x : op.GF_elem) z += K.Gcpt(x.c,x.r)*x.v;
    for(auto &x : op.IGF_elem) z += K.Gcpt(x.c,x.r)*x.v*K.phase[x.n] + K.Gcpt(x.r,x.c)*conjugate(x.v*K.phase[x.n]);
    I[0] = imag(z);
  };
  vector<double> Iv(1,0.0);
  QCM::k_integral(model->spatial_dimension, F, Iv, accur_OP, global_bool("verb_integrals"));

  double result = -M_1_PI*Iv[0]/model->Lc;

  if(model->mixing == HS_mixing::up_down){
    Green_function G_down = cluster_Green_function(w, false, true);
    G_down.w = w;
    G_down.spin_down = true;
    G_down.G.block.assign(n_clus, matrix<Complex>());
    for(size_t i = 0; i<n_clus; i++){
      if(model->clusters[i].conj){
        G_down.G.block[i] = matrix<Complex>(model->GF_dims[i], ED::Green_function(conjugate(w), true, n_clus*label+model->clusters[i].ref, false));
        G_down.G.block[i].cconjugate();
      }
      else
        G_down.G.block[i] = matrix<Complex>(model->GF_dims[i], ED::Green_function(w, true, n_clus*label+model->clusters[i].ref, false));
    }
    G_down.G.set_size();

    auto F_down = [this, G_down, op] (vector3D<double> &k, const int *nv, double *I) mutable {
      for(size_t i = 0 ; i<*nv; i++) I[i] = 0.0;
      Green_function_k K(G_down,k);
      set_Gcpt(K);
      Complex z(0.0);
      for(auto &x : op.GF_elem) z += K.Gcpt(x.c,x.r)*x.v;
      for(auto &x : op.IGF_elem) z += K.Gcpt(x.c,x.r)*x.v*K.phase[x.n] + K.Gcpt(x.r,x.c)*conjugate(x.v*K.phase[x.n]);
      I[0] = imag(z);
    };

    to_zero(Iv);
    QCM::k_integral(model->spatial_dimension, F_down, Iv, accur_OP, global_bool("verb_integrals"));
    result += -M_1_PI*Iv[0]/model->Lc;
  }
  if(model->mixing == HS_mixing::normal) result *= 2;

  return result;
}



//==============================================================================
vector<double> lattice_model_instance::momentum_profile(const lattice_operator& op, const vector<vector3D<double>> &k_set)
{
  double accur_OP = global_double("accur_OP");
  auto Fmp = [this, k_set, op] (Complex w, vector3D<double> &ki, const int *nv, double *I) mutable {
    const double w_offset = 2.0;
    Complex G_pole = 1.0/(w - w_offset);

    Green_function G = cluster_Green_function(w, false, false);
    size_t i = 0;
    for(auto& k : k_set){
      Green_function_k K(G,k);
      set_Gcpt(K);
      K.Gcpt.add(-G_pole); // regulates the Green function at high frequency (subtracts G_pole times the identity matrix)
      Complex z(0.0);
      for(auto &x : op.GF_elem) z += K.Gcpt(x.c,x.r)*x.v;
      for(auto &x : op.IGF_elem) z += K.Gcpt(x.c,x.r)*x.v*K.phase[x.n];
      I[i++] = real<double>(z);
    }
    if(model->mixing == HS_mixing::up_down){
      Green_function G_dw = cluster_Green_function(w,false,true);
      i = 0;
      for(auto& k : k_set){
        Green_function_k K(G_dw,k);
        set_Gcpt(K);
        K.Gcpt.add(-G_pole); // regulates the Green function at high frequency (subtracts G_pole times the identity matrix)
        Complex z(0.0);
        for(auto &x : op.GF_elem) z += K.Gcpt(x.c,x.r)*x.v;
        for(auto &x : op.IGF_elem) z += K.Gcpt(x.c,x.r)*x.v*K.phase[x.n];
        I[i++] += real<double>(z);
      }
    }
  };
  vector<double> A(k_set.size(), 0.0);
  QCM::wk_integral(0, Fmp, A, accur_OP, global_bool("verb_integrals"));
  if(model->mixing == HS_mixing::normal) A *= 2;
  else if((model->mixing&HS_mixing::anomalous)  == HS_mixing::anomalous) A += op.nambu_correction;
  else if(model->mixing == HS_mixing::full) {A += op.nambu_correction_full; A *= 0.5;}
  A *= op.norm;
  return A;
}

//==============================================================================
vector<double> lattice_model_instance::momentum_profile_per(const lattice_operator& op, const vector<vector3D<double>> &k_set)
{
  double accur_OP = global_double("accur_OP");
  auto Fmp = [this, k_set, op] (Complex w, vector3D<double> &ki, const int *nv, double *I) mutable {
    const double w_offset = 2.0;
    Complex G_pole = 1.0/(w - w_offset);

    Green_function G = cluster_Green_function(w, false, false);
    size_t i = 0;
    for(auto& k : k_set){
      vector3D<double> k_red = model->superdual.to(model->physdual.from(k));
      Green_function_k K(G,k_red);
      periodized_Green_function(K);
      K.g.add(-G_pole); // regulates the Green function at high frequency (subtracts G_pole times the identity matrix)
      matrix<Complex> S(G.G.r);
      // build the operator matrix
      for(auto &x : op.GF_elem) S(x.r, x.c) += x.v;
      for(auto &x : op.IGF_elem) S(x.r, x.c) += x.v*K.phase[x.n];
      matrix<Complex> S_per = model->periodize(k_red, S);
      I[i++] = realpart(K.g.trace_product(S_per));
    }
    if(model->mixing == HS_mixing::up_down){
      Green_function G_dw = cluster_Green_function(w,false,true);
      i = 0;
      for(auto& k : k_set){
        vector3D<double> k_red = model->superdual.to(model->physdual.from(k));
        Green_function_k K(G_dw,k_red);
        periodized_Green_function(K);
        K.g.add(-G_pole); // regulates the Green function at high frequency (subtracts G_pole times the identity matrix)
        matrix<Complex> S(G.G.r);
        // build the operator matrix
        for(auto &x : op.GF_elem) S(x.r, x.c) += x.v;
        for(auto &x : op.IGF_elem) S(x.r, x.c) += x.v*K.phase[x.n];
        matrix<Complex> S_per = model->periodize(k_red, S);
        I[i++] = realpart(K.g.trace_product(S_per));
      }
    }
  };
  vector<double> A(k_set.size(), 0.0);
  QCM::wk_integral(0, Fmp, A, accur_OP, global_bool("verb_integrals"));
  if(model->mixing == HS_mixing::normal) A *= 2;
  else if((model->mixing&HS_mixing::anomalous)  == HS_mixing::anomalous) A += op.nambu_correction;
  else if(model->mixing == HS_mixing::full) {A += op.nambu_correction_full; A *= 0.5;}
  A *= op.norm;
  return A;
}


#define WINF 1e5
//==============================================================================
double lattice_model_instance::potential_energy()
{
  if(PE_solved) return E_pot;
  PE_solved = true;

  double accur_OP = global_double("accur_OP");
  if(!gf_solved) Green_function_solve();

  // computing the infinite frequency limit
  double prev_cutoff = global_double("cutoff_scale");
  set_global_double("cutoff_scale", WINF); // important for convergence


  tr_sigma_inf = 0.0;
  for(int i=0; i<model->clusters.size(); i++){
    tr_sigma_inf += ED::tr_sigma_inf(i + label*model->clusters.size());
  }

  // frequency integral
  vector<double> I(1);
  I[0] = 0.0;

  if(grid_freqs.size() > 0){
    // ---- Discrete (omega, k) grid path ----
    // Mirrors the averages()/Potthoff_functional() pattern: build the cluster
    // Green function (and its self-energy) once per frequency, then parallelize
    // the k-sum via k_integral_grid. The cluster GF is captured by value into
    // the lambda and read concurrently inside the OpenMP region (no shared
    // mutable state).
    int nkx, nky, nkz;
    QCM::get_wavevector_grid(nkx, nky, nkz);
    if(global_bool("verb_integrals"))
      cout << "computing potential energy from a grid of " << grid_freqs.size() << " frequencies and "
           << nkx << "x" << nky << "x" << nkz << " k-points (dim=" << model->spatial_dimension << ")" << endl;
    bool has_dn = (model->mixing == HS_mixing::up_down);
    vector<double> Iw(1);
    for(int iw = 0; iw < (int)grid_freqs.size(); iw++){
      Complex wc(0, grid_freqs[iw]);
      Green_function G_up = cluster_Green_function(wc, true, false);
      Green_function G_dn;
      if(has_dn) G_dn = cluster_Green_function(wc, true, true);

      auto F_k = [this, G_up, G_dn, has_dn] (vector3D<double> &k, const int *nv, double *I) mutable {
        potential_energy_integrand_k(G_up, has_dn ? &G_dn : nullptr, k, nv, I);
      };
      Iw[0] = 0.0;
      QCM::k_integral_grid(model->spatial_dimension, nkx, nky, nkz, F_k, Iw);
      I[0] += Iw[0] * grid_weights[iw];
    }
  }
  else{
    // ---- Adaptive cubature in (omega, k) ----
    auto F = [this] (Complex w, vector3D<double> &k, const int *nv, double *I) mutable {potential_energy_integrand(w, k, nv, I);};
    QCM::wk_integral(model->spatial_dimension, F, I, 0.1*accur_OP, global_bool("verb_integrals"));
  }

  set_global_double("cutoff_scale", prev_cutoff); // restores to previous value
  if(model->mixing == HS_mixing::full) I[0] *= 0.5; // full Nambu doubling overestimates by a factor of 2
  E_pot = 0.5*I[0]/model->Lc;
  return E_pot;
}




//==============================================================================
complex<double> lattice_model_instance::TrSigmaG(Complex w, vector3D<double> &k, bool spin_down)
{
  complex<double> z(0.0);
  Green_function G = cluster_Green_function(w, false, spin_down);
  Green_function_k K(G,k);
  set_Gcpt(K);
  cluster_self_energy(G);
  z = G.sigma.build_matrix().trace_product(K.Gcpt);
  return z;
}

//==============================================================================
void lattice_model_instance::potential_energy_integrand(Complex w, vector3D<double> &k, const int *nv, double *I)
{
  check_signals();

  const double w_offset = 2.0;

  Complex e = TrSigmaG(w, k, false);
  if(model->mixing == HS_mixing::up_down) e += TrSigmaG(w, k, true);
  else if(model->mixing == HS_mixing::normal) e *= 2;

  e -= tr_sigma_inf/(w-w_offset);
  I[0] = real<double>(e);
}


//==============================================================================
void lattice_model_instance::potential_energy_integrand_k(Green_function &G_up, Green_function *G_down, vector3D<double> &k, const int *nv, double *I)
{
  const double w_offset = 2.0;
  Complex w = G_up.w;

  Green_function_k K(G_up, k);
  set_Gcpt(K);
  Complex e = G_up.sigma.build_matrix().trace_product(K.Gcpt);

  if(G_down){
    Green_function_k K_dn(*G_down, k);
    set_Gcpt(K_dn);
    e += G_down->sigma.build_matrix().trace_product(K_dn.Gcpt);
  }
  else if(model->mixing == HS_mixing::normal) e *= 2;

  e -= tr_sigma_inf/(w-w_offset);
  I[0] = real<double>(e);
}


//==============================================================================
matrix<complex<double>> lattice_model_instance::Green_integral(bool spin_down, int clus)
{
  size_t dim;
  if(clus) dim = model->GF_dims[clus-1];
  else dim = model->dim_GF;
  double accur_OP = global_double("accur_OP");
  auto Fmp = [this, dim, spin_down] (Complex w, vector3D<double> &ki, const int *nv, double *I) mutable {
    const double w_offset = 2.0;
    Complex G_pole = 1.0/(w - w_offset);
    G_pole += conjugate(G_pole);
      matrix<Complex> Gloc = projected_Green_function(w, spin_down);
    Gloc.symmetrize();
    Gloc.v *= 2.0;
    Gloc.add(-G_pole);
    for(size_t i = 0 ; i<*nv; i++) I[i] = 0.0;
    hermitian_matrix_to_real_vector(Gloc, I);
  };
  auto Fmp_clus = [this, dim, spin_down] (Complex w, vector3D<double> &ki, const int *nv, double *I) mutable {
    const double w_offset = 2.0;
    Complex G_pole = 1.0/(w - w_offset);
    G_pole += conjugate(G_pole);
      matrix<Complex> Gloc = cluster_Green_function(0, w, spin_down, false);
    Gloc.symmetrize();
    Gloc.v *= 2.0;
    Gloc.add(-G_pole);
    for(size_t i = 0 ; i<*nv; i++) I[i] = 0.0;
    hermitian_matrix_to_real_vector(Gloc, I);
  };

  vector<double> A(dim*dim, 0.0);
  if(clus) QCM::wk_integral(0, Fmp_clus, A, accur_OP, global_bool("verb_integrals"));
  else QCM::wk_integral(0, Fmp, A, accur_OP, global_bool("verb_integrals"));
  if(model->mixing == HS_mixing::normal) A *= 2;
  else if(model->mixing == HS_mixing::full) {A *= 0.5;}
  return hermitian_matrix_from_real_vector(dim, A);
}