obe_tb.hpp

Go to the documentation of this file.

// Copyright (c) 2025--present, The Simons Foundation
// This file is part of TRIQS/modest and is licensed under the terms of GPLv3 or later.
// SPDX-License-Identifier: GPL-3.0-or-later
// See LICENSE in the root of this distribution for details.
 
#pragma once
#include <mpi/mpi.hpp>
#include <nda/nda.hpp>
 
#include <triqs/experimental.hpp>
#include <triqs/mesh/imfreq.hpp>
#include <triqs/gfs.hpp>
#include <triqs/experimental/utility/root_finder.hpp>
#include <triqs/utility/macros.hpp>
 
#include "downfolding.hpp"
#include "loaders.hpp"
#include "triqs_modest/utils/gf_supp.hpp"
 
namespace triqs::modest {
 
  using namespace triqs::experimental::lattice;
  using namespace triqs::experimental::utility;

  struct one_body_elements_tb {
    local_space C_space;  
    std::vector<tb_hk> H; 
    // downfolding_projector P;
    // C2PY_IGNORE std::optional<ibz_symmetry_ops> ibz_symm_ops = {}; //< IBZ symmetrizer after a k-sum

    one_body_elements_tb(std::vector<tb_hk> H_sigma, local_space ls) : C_space{std::move(ls)}, H{std::move(H_sigma)} {};

    one_body_elements_tb(std::vector<tb_hk> H_sigma, spin_kind_e spin_kind, std::vector<atomic_orbs> atomic_shells);

    bool operator==(one_body_elements_tb const &) const = default;

    C2PY_IGNORE friend void mpi_broadcast(one_body_elements_tb &x, mpi::communicator c = {}, int root = 0) {
      mpi::broadcast(x.C_space, c, root);
      mpi::broadcast(x.H, c, root);
    }
  };


  one_body_elements_tb one_body_elements_from_wannier90(std::string const &wannier_file_path, spin_kind_e spin_kind,
                                                        std::vector<atomic_orbs> atomic_shells);

  one_body_elements_tb one_body_elements_from_wannier90(std::string const &wannier_file_path_up, std::string const &wannier_file_path_dn,
                                                        spin_kind_e spin_kind, std::vector<atomic_orbs> atomic_shells);
 // tb  factories

  nda::array<nda::matrix<dcomplex>, 2> Hloc(std::vector<tb_hk> const &H_sigma, std::vector<atomic_orbs> const &atomic_shells);

  nda::array<nda::matrix<dcomplex>, 2> impurity_levels(one_body_elements_tb const &obe);

  one_body_elements_tb fold(superlattice const &sl, one_body_elements_tb const &obe);

  one_body_elements_tb rotate(one_body_elements_tb const &obe, nda::matrix<dcomplex> const &U);
 
  one_body_elements_tb extend_to_spin(one_body_elements_tb const &obe);
  one_body_elements_tb add_local_term(one_body_elements_tb const &obe, nda::matrix<dcomplex> const &local_term);
 
  //  -----------------------------------------------------------------------


  template <typename Mesh>
  block2_gf<Mesh, matrix_valued> gloc(one_body_elements_tb const &obe, double mu, block2_gf<Mesh, matrix_valued> const &Sigma_dynamic,
                                      nda::array<nda::matrix<dcomplex>, 2> const &Sigma_static, bz_int_options const &opt) {
 
    auto &mesh       = Sigma_dynamic(0, 0).mesh();
    auto n_sigma     = obe.C_space.n_sigma();
    auto gloc_result = make_block2_gf(mesh, obe.C_space.Gc_block_shape());
    if (n_sigma != Sigma_static.shape(1)) { throw std::runtime_error("Mismatch between the spin channels in Sigma_Static and Sigma_Dynamic"); }
    if (n_sigma != obe.H.size()) { throw std::runtime_error("Mismatch between the spin channels in Sigma and spin channels in Hamiltonian."); }
 
    // Embedding decomposition from structure of Sigma -- provides a list of block names
    auto embedding_decomp = get_struct(Sigma_dynamic).dims(r_all, 0) | tl::to<std::vector>();
 
    for (auto sigma : range(n_sigma)) {
 
      // spin index
      auto Sigma_full_space = gfs::gf(mesh, {obe.H[sigma].n_orbitals(), obe.H[sigma].n_orbitals()}); //
      for (auto &&[block, R] : enumerated_sub_slices(embedding_decomp)) {
        for (auto [n, w] : enumerate(mesh)) {
          Sigma_full_space.data()(n, R, R) = Sigma_dynamic(block, sigma).data()(n, r_all, r_all) + Sigma_static(block, sigma);
        }
      }
      // Call the TRIQS version of this function
      gloc_result(0, sigma) = gloc(obe.H[sigma], mu, Sigma_full_space, opt);
    }
    return gloc_result;
  }

  template <typename Mesh>
  block2_gf<Mesh, matrix_valued> gloc(Mesh const &mesh, one_body_elements_tb const &obe, double mu, bz_int_options const &opt) {
    auto result = make_block2_gf(mesh, obe.C_space.Gc_block_shape());
    for (auto sigma : range(obe.C_space.n_sigma())) { result(0, sigma) = gloc(mesh, obe.H[sigma], mu, opt); }
    return result;
  }

 
  //  -----------------------------------------------------------------------

  template <typename Mesh>
  double density(one_body_elements_tb const &obe, double mu, block2_gf<Mesh, matrix_valued> const &Sigma_dynamic,
                 nda::array<nda::matrix<dcomplex>, 2> const &Sigma_static, bz_int_options const &opt) {
 
    // auto n_blocks = Sigma_dynamic.size1();
    auto n_sigma = obe.C_space.n_sigma();
 
    double n  = 0;
    auto Gloc = gloc(obe, mu, Sigma_dynamic, Sigma_static, opt); // returns block2gf (1, nsigma, {norb, norb})
    // return type of Gloc is a B2GF with dimensions (1, nspin, {norb, norb})
    for (auto sigma : range(n_sigma)) { // spin index
      n += real(nda::trace(density(Gloc(0, sigma))));
    }
    return n;
  }

  template <typename Mesh> double density(one_body_elements_tb const &obe, double mu, Mesh const &mesh, bz_int_options const &opt) {
 
    auto Sigma_dynamic = make_block2_gf(mesh, obe.C_space.Gc_block_shape());
    auto Sigma_static  = nda::array<nda::matrix<dcomplex>, 2>(1, obe.C_space.n_sigma());
    for (auto [i, j] : Sigma_static.indices()) { Sigma_static(i, j) = nda::zeros<dcomplex>(obe.C_space.dim(), obe.C_space.dim()); }
    return density(obe, mu, Sigma_dynamic, Sigma_static, opt);
  }
 
  //  -----------------------------------------------------------------------
  template <typename Mesh>
  double find_chemical_potential(double const target_density, one_body_elements_tb const &obe, block2_gf<Mesh, matrix_valued> const &Sigma_dynamic,
                                 nda::array<nda::matrix<dcomplex>, 2> const &Sigma_static, bz_int_options const &opt,
                                 std::string method = "dichotomy", double precision = 1.e-5, bool verbosity = true) {
    std::function<double(double)> f = [&obe, &Sigma_dynamic, &Sigma_static, &opt](double x) {
      return density(obe, x, Sigma_dynamic, Sigma_static, opt);
    };
    return std::get<0>(
       root_finder(method, f, 0.0, target_density, precision, 0.5, 1000, "Chemical Potential", "Total Density", verbosity));
  }

  template <typename Mesh>
  double find_chemical_potential(double const target_density, one_body_elements_tb const &obe, Mesh const &mesh, bz_int_options const &opt,
                                 std::string method = "dichotomy", double precision = 1.e-5, bool verbosity = true) {
 
    auto Sigma_dynamic = make_block2_gf(mesh, obe.C_space.Gc_block_shape());
    auto Sigma_static  = nda::array<nda::matrix<dcomplex>, 2>(1, obe.C_space.n_sigma());
    for (auto [i, j] : Sigma_static.indices()) { Sigma_static(i, j) = nda::zeros<dcomplex>(obe.C_space.dim(), obe.C_space.dim()); }
    return find_chemical_potential(target_density, obe, Sigma_dynamic, Sigma_static, opt, method, precision, verbosity);
  }
 
  // --------  instantiations --------------

  template block2_gf<mesh::imfreq, matrix_valued> gloc(one_body_elements_tb const &obe, double mu,
                                                       block2_gf<mesh::imfreq, matrix_valued> const &Sigma_dynamic,
                                                       nda::array<nda::matrix<dcomplex>, 2> const &Sigma_static, bz_int_options const &opt);
 
  template block2_gf<mesh::imfreq, matrix_valued> gloc(mesh::imfreq const &mesh, one_body_elements_tb const &obe, double mu,
                                                       bz_int_options const &opt);
 
  template double density(one_body_elements_tb const &obe, double mu, block2_gf<mesh::imfreq, matrix_valued> const &Sigma_dynamic,
                          nda::array<nda::matrix<dcomplex>, 2> const &Sigma_static, bz_int_options const &opt);
 
  template double density(one_body_elements_tb const &obe, double mu, mesh::imfreq const &mesh, bz_int_options const &opt);
 
  template double find_chemical_potential(double const target_density, one_body_elements_tb const &obe,
                                          block2_gf<mesh::imfreq, matrix_valued> const &Sigma_dynamic,
                                          nda::array<nda::matrix<dcomplex>, 2> const &Sigma_static, bz_int_options const &opt, std::string method,
                                          double precision, bool verbosity);
 
  template double find_chemical_potential(double const target_density, one_body_elements_tb const &obe, mesh::imfreq const &mesh,
                                          bz_int_options const &opt, std::string method, double precision, bool verbosity);

 
} // namespace triqs::modest