Program Listing for File QCM.cpp

↰ Return to documentation for file (src_qcm/QCM.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;
  }

}