.. _program_listing_file_src_qcm_average.cpp:

Program Listing for File average.cpp
====================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_qcm_average.cpp>` (``src_qcm/average.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: 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);
   }