.. _program_listing_file_src_qcm_QCM.cpp:

Program Listing for File QCM.cpp
================================

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

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

.. code-block:: cpp

   //
   //  QCM.cpp
   //  qcm3
   //
   //  Created by David Sénéchal on 16-11-22.
   //  Copyright © 2016 David Sénéchal. All rights reserved.
   //
   
   #include <memory>
   #include <iostream>
   #include <fstream>
   #ifdef _OPENMP
     #include <omp.h>
   #endif
   
   #include "QCM.hpp"
   #include "global_parameter.hpp"
   #include "lattice_model.hpp"
   #include "lattice_model_instance.hpp"
   #include "parser.hpp"
   #include "qcm_ED.hpp"
   #include "parameter_set.hpp"
   #include "model.hpp"
   #include "model_instance_base.hpp"
   
   //==============================================================================
   // global variables
   
   shared_ptr<lattice_model> qcm_model = nullptr; // pointer to the unique lattice model of the library
   map<int, unique_ptr<lattice_model_instance>> lattice_model_instances; // list of instances
   
   extern map<string, shared_ptr<model>> models;
   extern map<size_t, shared_ptr<model_instance_base>> model_instances;
   
   vector<double> grid_freqs; // optional imaginary frequency grid for discrete integrals
   vector<double> grid_weights; // weights associated with grid_freqs
   
   // per-direction wavevector grid sizes set via QCM::set_wavevector_grid().
   // A value of 0 means "not set: fall back to global_int(\"kgrid_side\")".
   int wavevector_grid_nx = 0;
   int wavevector_grid_ny = 0;
   int wavevector_grid_nz = 0;
   
   //==============================================================================
   namespace QCM{
   
   
   void set_wavevector_grid(int nkx, int nky, int nkz){
     wavevector_grid_nx = nkx;
     wavevector_grid_ny = nky;
     wavevector_grid_nz = nkz;
   }
   
   
   void get_wavevector_grid(int& nkx, int& nky, int& nkz){
     int def = (int)global_int("kgrid_side");
     nkx = wavevector_grid_nx > 0 ? wavevector_grid_nx : def;
     nky = wavevector_grid_ny > 0 ? wavevector_grid_ny : def;
     nkz = wavevector_grid_nz > 0 ? wavevector_grid_nz : def;
   }
   
   
   void great_reset(){
     for(auto& x:lattice_model_instances) x.second.reset();
     lattice_model_instances.clear();
     for(auto& x:models) x.second.reset();
     models.clear();
     qcm_model->sector_strings.clear();
     qcm_model->param_set.reset();
     parameter_set::parameter_set_defined = false;
     lattice_model::model_consolidated = false;
     qcm_model.reset();
     qcm_model = make_shared<lattice_model>();
   }
   
   
   void erase_lattice_model_instance(size_t label){
       size_t nsys = qcm_model->nsys;
       for(size_t i=0; i < nsys; i++){
         model_instances.erase(label*nsys+i);
       }
       lattice_model_instances.erase(label);
   }
   
   
   
   
     void check_instance(int label)
     {
       if(lattice_model_instances.find(label) == lattice_model_instances.end()) 
         qcm_throw("The instance # "+to_string(label)+" does not exist.");
     }
   
   
   
     void print_model(const string& filename, bool asy_operators, bool asy_labels, bool asy_orb, bool asy_neighbors, bool asy_working_basis)
     {
       ofstream fout(filename);
       if (!fout.good()) qcm_throw("failed to open file " + filename);
       qcm_model->print(fout, asy_operators, asy_labels, asy_orb, asy_neighbors, asy_working_basis);
       fout.close();
     }
     
     
     
     
     void new_model_instance(int label)
     {
       if(qcm_model->param_set == nullptr) qcm_throw("The parameters have not been specified yet.");
       lattice_model_instances[label] = unique_ptr<lattice_model_instance>(new lattice_model_instance(qcm_model, label));
     }
     
   
     
     
     void set_parameters(vector<pair<string,double>>& val, vector<tuple<string, double, string>>& equiv)
     {
       qcm_model->param_set = make_shared<parameter_set>(qcm_model, val, equiv);
     }
   
     
     void set_parameter(const string& name, double value)
     {
       if(qcm_model->param_set == nullptr) qcm_throw("The parameters have not been specified yet.");
       qcm_model->param_set->set_value(name, value);
     }
     
     
     
     void set_multiplier(const string& name, double value)
     {
       if(qcm_model->param_set == nullptr) qcm_throw("The parameters have not been specified yet.");
       qcm_model->param_set->set_multiplier(name, value);
     }
     
     
     
     map<string,double> parameters()
     {
       if(qcm_model->param_set == nullptr) qcm_throw("The parameters have not been specified yet.");
       return qcm_model->param_set->value_map();
     }
     
     
     map<string,double> instance_parameters(int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->params;
     }
   
     
     
     vector<double> momentum_profile(const string& op, const vector<vector3D<double>> &k_set, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       if(qcm_model->term.find(op) == qcm_model->term.end()) qcm_throw("The lattice operator "+op+" does not exist.");
       return lattice_model_instances.at(label)->momentum_profile_per(*qcm_model->term.at(op), k_set);
     }
     
     
     
     vector<pair<double,string>> ground_state(int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->ground_state();
     }
     
     
     
     matrix<complex<double>> cluster_Green_function(size_t i, complex<double> w, bool spin_down, int label, bool blocks)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->cluster_Green_function(i, w, spin_down, blocks);
     }
     
     
     matrix<complex<double>> cluster_self_energy(size_t i, complex<double> w, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->cluster_self_energy(i, w, spin_down);
     }
     
     
   
     matrix<complex<double>> cluster_hopping_matrix(size_t i, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->cluster_hopping_matrix(i, spin_down);
     }
     
     
     
     
     
     matrix<complex<double>> hybridization_function(complex<double> w, bool spin_down, size_t i, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->hybridization_function(i, w, spin_down);
     }
     
     
   
   
     matrix<complex<double>> Green_integral(bool spin_down, int clus, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->Green_integral(spin_down, clus);
     }
     
   
   
   
   
   
     matrix<complex<double>> CPT_Green_function(const complex<double> w, const vector3D<double> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances.at(label)->model;
       vector3D<double> K = mod.superdual.to(mod.physdual.from(k));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       Green_function_k M(G, K);
       lattice_model_instances.at(label)->set_Gcpt(M);
       return M.Gcpt;
     }
     
     
     vector<matrix<complex<double>>> CPT_Green_function(const complex<double> w, const vector<vector3D<double>> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances[label]->model;
       vector<vector3D<double>> K(k.size());
       for(size_t i = 0; i< K.size(); i++) K[i] = mod.superdual.to(mod.physdual.from(k[i]));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       vector<matrix<Complex>> R(k.size());
       for(size_t i = 0; i< k.size(); i++) {
         Green_function_k M(G, K[i]);
         lattice_model_instances.at(label)->set_Gcpt(M);
         R[i] = M.Gcpt;
       }
       return R;
     }
     
     
     matrix<complex<double>> CPT_Green_function(int iw, int ik, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances.at(label)->model;
       if(mod.hybrid == nullptr) qcm_throw("CPT_Green_function(iw, ik, ...) requires an external hybridization (hybrid_file)");
       lattice_hybrid& H = *mod.hybrid;
       if(iw < 0 || (size_t)iw >= H.nw) qcm_throw("frequency index iw out of range");
       if(ik < 0 || (size_t)ik >= H.nk) qcm_throw("wavevector index ik out of range");
       Complex w(0.0, H.w[iw]);
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, false);
       G.iw = iw;
       Green_function_k M(G, H.k[ik], ik);
       lattice_model_instances.at(label)->set_Gcpt(M);
       cout << "k (read) = " << H.k[ik] << "\tw = " << w << endl; // TEMPO
       cout << "gamma = \n" << mod.lattice_hybridization(M.G.iw, M.ik) << endl; // TEMPO
       return M.Gcpt;
     }
   
   
     matrix<complex<double>> V_matrix(const complex<double> w, const vector3D<double> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances.at(label)->model;
       vector3D<double> K = mod.superdual.to(mod.physdual.from(k));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       Green_function_k M(G, K);
       lattice_model_instances.at(label)->set_V(M);
       return M.V;
     }
     
     
     vector<matrix<complex<double>>> CPT_Green_function_inverse(const complex<double> w, const vector<vector3D<double>> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances[label]->model;
       vector<vector3D<double>> K(k.size());
       for(size_t i = 0; i< K.size(); i++) K[i] = mod.superdual.to(mod.physdual.from(k[i]));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       G.G.inverse();
       vector<matrix<Complex>> R(k.size());
       for(size_t i = 0; i< k.size(); i++) {
         Green_function_k M(G, K[i]);
         lattice_model_instances.at(label)->inverse_Gcpt(G.G, M);
         R[i] = M.V;
       }
       return R;
     }
   
   
     matrix<complex<double>> tk(const vector3D<double> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances.at(label)->model;
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(Complex(0,0.1), false, spin_down);
       Green_function_k M(G, k);
       lattice_model_instances.at(label)->set_V(M);
       return M.t;
     }
     
     
     vector<matrix<complex<double>>> tk(const vector<vector3D<double>> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances[label]->model;
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(Complex(0., 0.1), false, spin_down);
       vector<matrix<Complex>> R(k.size());
       for(size_t i = 0; i< k.size(); i++) {
         Green_function_k M(G, k[i]);
         lattice_model_instances.at(label)->set_Gcpt(M);
         R[i] = M.t;
       }
       return R;
     }
   
     
     vector<double> dos(const complex<double> w, int label, bool use_grid)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->dos(w, use_grid);
     }
     
     
     double spectral_average(const string& name, const complex<double> w, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->spectral_average(name, w);
     }
   
     
     
     matrix<complex<double>> periodized_Green_function(const complex<double> w, const vector3D<double> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       vector3D<double> K = qcm_model->superdual.to(qcm_model->physdual.from(k));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       Green_function_k M(G, K);
       lattice_model_instances.at(label)->periodized_Green_function(M);
       return M.g;
     }
     
     
     matrix<complex<double>> band_Green_function(const complex<double> w, const vector3D<double> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       vector3D<double> K = qcm_model->superdual.to(qcm_model->physdual.from(k));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       Green_function_k M(G, K);
       lattice_model_instances.at(label)->periodized_Green_function(M);
       
       return lattice_model_instances.at(label)->band_Green_function(M);
     }
     
     
     
     
     vector<matrix<complex<double>>> periodized_Green_function(const complex<double> w, const vector<vector3D<double>> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       vector<vector3D<double>> K(k.size());
       for(size_t i = 0; i< K.size(); i++) K[i] = qcm_model->superdual.to(qcm_model->physdual.from(k[i]));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       vector<matrix<Complex>> R;
       R.reserve(K.size());
       for(size_t i = 0; i< K.size(); i++) {
         Green_function_k M(G, K[i]);
         lattice_model_instances.at(label)->periodized_Green_function(M);
         R.push_back(M.g);
       }
       return R;
     }
     
   
     vector<matrix<complex<double>>> band_Green_function(const complex<double> w, const vector<vector3D<double>> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       vector<vector3D<double>> K(k.size());
       for(size_t i = 0; i< K.size(); i++) K[i] = qcm_model->superdual.to(qcm_model->physdual.from(k[i]));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       vector<matrix<Complex>> R;
       R.reserve(K.size());
       for(size_t i = 0; i< K.size(); i++) {
         Green_function_k M(G, K[i]);
         lattice_model_instances.at(label)->periodized_Green_function(M);
         R.push_back(lattice_model_instances.at(label)->band_Green_function(M));
       }
       return R;
     }
     
     
   vector<complex<double>> periodized_Green_function_element(int r, int c, const complex<double> w, const vector<vector3D<double>> &k, bool spin_down, int label)
   {
     #ifdef QCM_DEBUG
   check_instance(label);
   #endif
     vector<vector3D<double>> K(k.size());
     for(size_t i = 0; i< K.size(); i++) K[i] = qcm_model->superdual.to(qcm_model->physdual.from(k[i]));
     Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
     vector<Complex> R;
     R.reserve(K.size());
     for(size_t i = 0; i< K.size(); i++) {
       Green_function_k M(G, K[i]);
       lattice_model_instances.at(label)->periodized_Green_function(M);
       R.push_back(M.g(r,c));
     }
     return R;
   }
     
     
     
   
     vector<vector<double>> dispersion(const vector<vector3D<double>> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances.at(label)->model;
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model_instance& inst = *lattice_model_instances.at(label);
       
       vector<vector3D<double>> K(k.size());
       for(size_t i = 0; i< K.size(); i++) K[i] = mod.superdual.to(mod.physdual.from(k[i]));
       Green_function G = inst.cluster_Green_function({0.0,0.05}, false, spin_down);
       
       vector<vector<double>> R(K.size());
       for(size_t i = 0; i< K.size(); i++) {
         Green_function_k M(G, K[i]);
         R[i] = inst.dispersion(M);
       }
       return R;
     }
     
     
     
     
     matrix<complex<double>> compact_tiling(const matrix<complex<double>>& A, const vector3D<double>& k)
     {
       lattice_model& mod = *qcm_model;
       return mod.compact_tiling(A, k);
     }
   
   
     ED::CombinedMCF_data get_combined_mcf_k(const vector3D<double>& k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances.at(label)->model;
       QCM_ASSERT(mod.clusters.size() == 1);
   
       // Convert to superdual basis (same convention as set_V / periodize)
       vector3D<double> K = mod.superdual.to(mod.physdual.from(k));
   
       // Build a dummy Green_function_k to retrieve V(k)
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(Complex(0, 0.1), false, spin_down);
       Green_function_k M(G, K);
       lattice_model_instances.at(label)->set_V(M);
   
       // Get combined MCF from the ED layer
       size_t sys_label = mod.nsys * label + mod.clusters[0].sys_start;
       auto D = ED::get_combined_mcf(spin_down, sys_label);
   
       // Add inter-cluster hopping into the first diagonal block
       D.A[0].v += M.V.v;
   
       int nA = (int)D.A.size();
       int nB = (int)D.B.size();
   
       // Apply compact_tiling to A[j>=1] and B[j>=1]
       for(int j = 1; j < nA; ++j)
         D.A[j] = mod.compact_tiling(D.A[j], K);
       for(int j = 1; j < nB; ++j)
         D.B[j] = mod.compact_tiling(D.B[j], K);
   
       // Periodize all blocks into the band basis
       ED::CombinedMCF_data Dper;
       Dper.A.resize(nA);
       Dper.B.resize(nB);
       for(int j = 0; j < nA; ++j)
         Dper.A[j] = mod.periodize(K, D.A[j]);
       for(int j = 0; j < nB; ++j)
         Dper.B[j] = mod.periodize(K, D.B[j]);
       Dper.W = mod.periodize(K, D.W);
   
       return Dper;
     }
   
   
     vector<matrix<Complex>> epsilon(const vector<vector3D<double>> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances.at(label)->model;
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model_instance& inst = *lattice_model_instances.at(label);
       
       vector<vector3D<double>> K(k.size());
       for(size_t i = 0; i< K.size(); i++) K[i] = mod.superdual.to(mod.physdual.from(k[i]));
       Green_function G = inst.cluster_Green_function({0.0,0.05}, false, spin_down);
       
       vector<matrix<Complex>> R(K.size());
       for(size_t i = 0; i< K.size(); i++) {
         Green_function_k M(G, K[i]);
         R[i] = inst.epsilon(M);
       }
       return R;
     }
   
   
     
     vector<matrix<complex<double>>> self_energy(const complex<double> w, const vector<vector3D<double>> &k, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model& mod = *lattice_model_instances.at(label)->model;
       vector<vector3D<double>> K(k.size());
       for(size_t i = 0; i< K.size(); i++) K[i] = mod.superdual.to(mod.physdual.from(k[i]));
       Green_function G = lattice_model_instances.at(label)->cluster_Green_function(w, false, spin_down);
       vector<matrix<Complex>> R;
       R.reserve(K.size());
       for(size_t i = 0; i< K.size(); i++) {
         Green_function_k M(G, K[i]);
         lattice_model_instances.at(label)->self_energy(M);
         R.push_back(M.sigma);
       }
       return R;
     }
     
     
     
     
     double Potthoff_functional(int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->Potthoff_functional();
     }
     
     
     
     double potential_energy(int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->potential_energy();
     }
     
     
     
     double kinetic_energy(int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->E_kin;
     }
   
   
   
     vector<pair<string,double>> averages(const vector<string> &_ops, int label)
     {
       return lattice_model_instances.at(label)->averages(_ops);
     }
     
     
     
     matrix<complex<double>> projected_Green_function(const complex<double> w, bool spin_down, int label)
     {
       return lattice_model_instances.at(label)->projected_Green_function(w, spin_down);
     }
     
     
     
     size_t Green_function_dimension()
     {
       return qcm_model->dim_GF;
     }
     
     
     size_t reduced_Green_function_dimension()
     {
       char periodization = global_char("periodization");
       if(periodization == 'N') return qcm_model->dim_GF;
       else return qcm_model->dim_reduced_GF;
     }
     
     
     int spatial_dimension()
     {
       return qcm_model->spatial_dimension;
     }
     
     
     
     int mixing()
     {
       return qcm_model->mixing;
     }
     
     
   
     double Berry_flux(vector<vector3D<double>>& k, int orb, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->Berry_flux(k, orb, false);
     }
   
   
   
   
   
     vector<double> Berry_curvature(vector3D<double>& k1, vector3D<double>& k2, int nk, int orb, bool rec, int dir, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->Berry_curvature(k1, k2, nk, orb, rec, dir);
     }
     
     
     
     void add_cluster(const vector3D<int64_t> &cpos, const vector<vector3D<int64_t>> &pos, int ref, bool conj)
     {
       if(qcm_model->is_closed){
         qcm_warning("model already created. Ignoring.");
         return;
       }
       for(size_t i=0; i<pos.size(); i++) qcm_model->sites.push_back({qcm_model->clusters.size(), i, 0, pos[i]+cpos});
       if(ref == -1) ref = qcm_model->clusters.size();
       qcm_model->clusters.push_back({pos.size(), qcm_model->sites.size(), cpos, ref, 0, conj, 0, 0});
       // n_sites, n_bath, offset, name, position, ref, mixing, sys_start, nsys
     }
     
     
     void add_system(const string &name, const int clus)
     {
       if(qcm_model->is_closed){
         qcm_warning("model already created. Ignoring.");
         return;
       }
       if(clus > qcm_model->clusters.size())
         qcm_throw("This system cannot be added to the cluster (out of range, or maybe clusters not yet added to lattice model)");
     
       auto tmp = ED::model_size(name);
       if(get<0>(tmp) != qcm_model->clusters[clus].n_sites) qcm_throw("The number of sites of cluster "+name+" is inconsistent with the cluster model");
       if(get<1>(tmp) > 0) qcm_model->bath_exists = true;
       // n_sites, n_bath, name, clus, mixing, n_sym
       qcm_model->systems.push_back({get<0>(tmp), get<1>(tmp), name, clus, 0, get<2>(tmp)});
     }
     
     
   
     void new_lattice_model(const string &name, vector<int64_t> &superlattice, vector<int64_t> &unit_cell, const string &latt_hybrid)
     {
       if(qcm_model==nullptr) qcm_throw("no cluster has been added to the model!");
       if(qcm_model->is_closed){
         qcm_warning("model already created. Ignoring.");
         return;
       }
       qcm_model->name = name;
       qcm_model->superlattice = lattice3D(superlattice);
       qcm_model->unit_cell = lattice3D(unit_cell);
       qcm_model->phys.trivial();
       qcm_model->hybrid_file = latt_hybrid;
       qcm_model->hybrid = nullptr;
       qcm_model->pre_operator_consolidate();
     }
     
     
     void set_basis(vector<double> &basis)
     {
       qcm_model->phys = basis3D(basis);
       qcm_model->phys.inverse();
       qcm_model->phys.dual(qcm_model->physdual);
       qcm_model->physdual.init();
     }
     
     
     void interaction_operator(const string &name, vector3D<int64_t> &link, double amplitude, int orb1, int orb2, const string &type)
     {
       if(qcm_model->is_closed){
         qcm_warning("model already created and closed. Ignoring operator creation.");
         return;
       }
       qcm_model->interaction_operator(name, link, amplitude, orb1-1, orb2-1, type);
     }
     
     
     void hopping_operator(const string &name, vector3D<int64_t> &link, double amplitude, int orb1, int orb2, int tau, int sigma)
     {
       if(qcm_model->is_closed){
         qcm_warning("model already created and closed. Ignoring operator creation.");
         return;
       }
       qcm_model->hopping_operator(name, link, amplitude, orb1-1, orb2-1, tau, sigma);
     }
     
     
     void current_operator(const string &name, vector3D<int64_t> &link, double amplitude, int orb1, int orb2, int dir, bool pau)
     {
       if(qcm_model->is_closed){
         qcm_warning("model already created and closed. Ignoring operator creation.");
         return;
       }
       qcm_model->current_operator(name, link, amplitude, orb1-1, orb2-1, dir, pau);
     }
   
   
     void anomalous_operator(const string &name, vector3D<int64_t> &link, complex<double> amplitude, int orb1, int orb2, const string& type)
     {
       if(qcm_model->is_closed){
         qcm_warning("model already created and closed. Ignoring operator creation.");
         return;
       }
       qcm_model->anomalous_operator(name, link, amplitude, orb1-1, orb2-1, type);
     }
     
     void density_wave(const string &name, vector3D<int64_t> &link, complex<double> amplitude, int orb, vector3D<double> Q, double phase, const string& type)
     {
       if(qcm_model->is_closed){
         qcm_warning("model already created and closed. Ignoring operator creation.");
         return;
       }
       qcm_model->density_wave(name, link, amplitude, orb-1, Q, phase, type);
     }
   
     void explicit_operator(const string &name, const string &type, const vector<tuple<vector3D<int64_t>, vector3D<int64_t>, complex<double>>> &elem, int tau, int sigma)
     {
       if(qcm_model->is_closed){
         qcm_warning("model already created and closed. Ignoring operator creation.");
         return;
       }
       qcm_model->explicit_operator(name, type, elem, tau, sigma);
     }
   
     
     vector<tuple<string, int, int, int, int>> cluster_info()
     {
       vector<tuple<string, int, int, int, int>> info(qcm_model->clusters.size());
       for(int i=0; i<info.size(); i++){
         get<0>(info[i]) = (int)qcm_model->clusters[i].n_sites;
         get<1>(info[i]) = (int)qcm_model->clusters[i].nsys;
         get<2>(info[i]) = (int)qcm_model->clusters[i].sys_start;
         get<3>(info[i]) = (int)qcm_model->clusters[i].n_sites*qcm_model->n_mixed;
       }
       return info;
     }
   
   
   
     vector<tuple<string, int, int, int>> systems_info()
     {
       vector<tuple<string, int, int, int>> info(qcm_model->systems.size());
       for(int i=0; i<info.size(); i++){
         get<0>(info[i]) = qcm_model->systems[i].name;
         get<1>(info[i]) = (int)qcm_model->systems[i].n_sites;
         get<2>(info[i]) = (int)qcm_model->systems[i].n_bath;
         get<3>(info[i]) = (int)qcm_model->systems[i].clus;
       }
       return info;
     }
   
   
   
     pair<vector<array<double,9>>, vector<array<complex<double>, 11>>> site_and_bond_profile(int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       auto& M = *lattice_model_instances[label];
       return M.site_and_bond_profile();
     }
   
     
   
     vector<pair<vector<double>, vector<double>>> Lehmann_Green_function(vector<vector3D<double>> &k, int orb, bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       for(size_t i = 0; i< k.size(); i++) k[i] = qcm_model->superdual.to(qcm_model->physdual.from(k[i]));
       return lattice_model_instances[label]->Lehmann_Green_function(k, orb, spin_down);
     }
   
   
     void Green_function_solve(int label){
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model_instances.at(label)->Green_function_solve();  
     }
   
     void CDMFT_variational_set(vector<vector<string>>& varia){
       if(qcm_model->param_set == nullptr) qcm_throw("The parameters have not been specified yet.");
       qcm_model->param_set->CDMFT_variational_set(varia);
     }
   
     void CDMFT_host(const vector<double>& freqs, const vector<double>& weights, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model_instances.at(label)->CDMFT_host(freqs, weights); 
     }
   
     void set_CDMFT_host(int label, const vector<double>& freqs, const int clus, const vector<matrix<Complex>>& H, const bool spin_down)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       lattice_model_instances.at(label)->set_CDMFT_host(freqs, clus, H, spin_down);
     }
   
   
     vector<vector<matrix<Complex>>> get_CDMFT_host(bool spin_down, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->get_CDMFT_host(spin_down); 
     }
   
     double CDMFT_distance(const vector<double>& p, int clus, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->CDMFT_distance(p, clus);
     }
   
     vector<double> CDMFT_residuals(const vector<double>& p, int clus, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->CDMFT_residuals(p, clus);
     }
   
     vector<double> CDMFT_gradient(const vector<double>& p, int clus, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->CDMFT_gradient(p, clus);
     }
   
     double monopole(vector3D<double>& k, double a, int nk, int orb, bool rec, int label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->monopole(k, a, nk, orb, rec); 
     }
   
   
     bool complex_HS(size_t label)
     {
       #ifdef QCM_DEBUG
       check_instance(label);
       #endif
       return lattice_model_instances.at(label)->complex_HS;
     }
   
   }