tight_binding.hpp

Go to the documentation of this file.

// Copyright (c) 2013-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2013-2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2016 Igor Krivenko
// Copyright (c) 2018-2023 Simons Foundation
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You may obtain a copy of the License at
//     https://www.gnu.org/licenses/gpl-3.0.txt
//
// Authors: Michel Ferrero, Igor Krivenko, Olivier Parcollet, Nils Wentzell, Thomas Ayral

 
#pragma once
 
#include "./brillouin_zone.hpp"
#include "../gfs/gf/gf.hpp"
#include "../mesh/brzone.hpp"
#include "../utility/macros.hpp"
 
#include <h5/h5.hpp>
#include <itertools/itertools.hpp>
#include <nda/linalg.hpp>
#include <nda/nda.hpp>
#include <nda/stdutil/complex.hpp>
 
#include <cmath>
#include <iostream>
#include <string>
#include <utility>
#include <vector>
 
namespace triqs::lattice {


  struct C2PY_IGNORE hopping_dict {
    std::vector<nda::vector<long>> displ_vec;

    std::vector<nda::matrix<dcomplex>> overlap_mat_vec;
  };

  class tight_binding {
 
    bravais_lattice bl_;
    std::vector<nda::vector<long>> displ_vec_;
    std::vector<nda::matrix<dcomplex>> overlap_mat_vec_;
 
    void check_hoppings();
 
    public:
    tight_binding(bravais_lattice bl, std::vector<nda::vector<long>> displ_vec, std::vector<nda::matrix<dcomplex>> overlap_mat_vec);

    tight_binding(bravais_lattice bl, hopping_dict hoppings);

    bravais_lattice const &lattice() const { return bl_; }

    auto const &displ_vec() const { return displ_vec_; }

    auto const &overlap_mat_vec() const { return overlap_mat_vec_; }

    template <typename R> r_t lattice_to_real_coordinates(R const &x) const { return bl_.lattice_to_real_coordinates(x); }

    long n_orbitals() const { return bl_.n_orbitals(); }

    template <typename F> friend void foreach (tight_binding const &tb, F f) {
      int n = static_cast<int>(tb.displ_vec_.size());
      for (int i = 0; i < n; ++i) f(tb.displ_vec_[i], tb.overlap_mat_vec_[i]);
    }

    template <typename K>
      requires(nda::ArrayOfRank<K, 1> or nda::ArrayOfRank<K, 2>)
    auto fourier(K const &k) const {
      // Make sure to account for ndim in lattice
      auto k_ndim = make_regular(k(nda::ellipsis(), range(lattice().ndim())));
 
      auto vals = [&](int j) {
        if constexpr (nda::ArrayOfRank<K, 1>) {
          return std::exp(2i * M_PI * nda::linalg::dot_generic(k_ndim, displ_vec_[j])) * overlap_mat_vec_[j];
        } else { // Rank==2
          auto k_mat = nda::make_matrix_view(k_ndim);
          auto exp   = [](auto d) { return std::exp(d); };
          auto exp_j = make_regular(nda::map(exp)(2i * M_PI * k_mat * displ_vec_[j]));
          return nda::linalg::outer_product(exp_j, overlap_mat_vec_[j]);
        }
      };
      auto res = make_regular(vals(0));
      for (int i = 1; i < displ_vec_.size(); ++i) res += vals(i);
      return res;
    }

    inline auto fourier(mesh::brzone const &k_mesh) const {
      auto kvecs = nda::matrix<double>(k_mesh.size(), 3);
      for (auto [n, k] : itertools::enumerate(k_mesh)) { kvecs(n, range::all) = k.value(); }
      auto kvecs_rec = make_regular(kvecs * k_mesh.bz().reciprocal_matrix_inv());
      auto h_k       = gfs::gf<mesh::brzone, gfs::matrix_valued>(k_mesh, {n_orbitals(), n_orbitals()});
      h_k.data()     = fourier(kvecs_rec);
      return h_k;
    }

    inline auto fourier(int n_l) const {
      auto k_mesh = mesh::brzone(brillouin_zone{bl_}, n_l);
      return fourier(k_mesh);
    }

    template <typename K>
    auto dispersion(K const &k) const
      requires(nda::ArrayOfRank<K, 1> or nda::ArrayOfRank<K, 2>)
    {
      if constexpr (nda::ArrayOfRank<K, 1>) {
        return nda::linalg::eigvalsh(fourier(k));
      } else { // Rank==2
        auto h_k = fourier(k);
        auto n_k = h_k.shape()[0];
        auto res = nda::array<double, 2>(n_k, n_orbitals());
        for (auto l : range(n_k)) res(l, range::all) = nda::linalg::eigvalsh(h_k(l, nda::ellipsis()));
        return res;
      }
    }

    inline auto dispersion(mesh::brzone const &k_mesh) const {
      auto h_k = fourier(k_mesh);
      auto e_k = gfs::gf<mesh::brzone, gfs::tensor_real_valued<1>>(k_mesh, {n_orbitals()});
      for (auto k : k_mesh) e_k[k] = nda::linalg::eigvalsh(h_k[k]);
      return e_k;
    }

    inline auto dispersion(int n_l) const {
      auto k_mesh = mesh::brzone(brillouin_zone{bl_}, n_l);
      return dispersion(k_mesh);
    }
 
    // ------------------- Comparison -------------------

    bool operator==(tight_binding const &tb) const { return bl_ == tb.bl_ && overlap_mat_vec_ == tb.overlap_mat_vec_ && displ_vec_ == tb.displ_vec_; }

    bool operator!=(tight_binding const &tb) const { return !(operator==(tb)); }
 
    // -------------------- print -------------------

    friend std::ostream &operator<<(std::ostream &sout, tight_binding const &tb) {
      sout << "Tight Binding Hamiltonian on " << tb.lattice() << "\nwith hoppings [";
      for (auto const &[displ, overlap_mat] : itertools::zip(tb.displ_vec(), tb.overlap_mat_vec())) sout << "\n   " << displ << " : " << overlap_mat;
      return sout << " ]";
    }
 
    // ------------------- HDF5 Read / Write -------------------

    [[nodiscard]] static std::string hdf5_format() { return "tight_binding"; }

    friend void h5_write(h5::group fg, std::string subgroup_name, tight_binding const &tb) {
      auto grp = fg.create_group(subgroup_name);
      write_hdf5_format(grp, tb); // NOLINT
      h5_write(grp, "bravais_lattice", tb.bl_);
      h5_write(grp, "displ_vec", tb.displ_vec_);
      h5_write(grp, "overlap_mat_vec", tb.overlap_mat_vec_);
    }

    CPP2PY_IGNORE
    static tight_binding h5_read_construct(h5::group g, std::string subgroup_name) {
      auto grp             = g.open_group(subgroup_name);
      auto bl              = h5::h5_read<bravais_lattice>(grp, "bravais_lattice");
      auto displ_vec       = h5::h5_read<std::vector<nda::vector<long>>>(grp, "displ_vec");
      auto overlap_mat_vec = h5::h5_read<std::vector<nda::matrix<dcomplex>>>(grp, "overlap_mat_vec");
      return tight_binding(bl, displ_vec, overlap_mat_vec);
    }
 
  }; // tight_binding

  std::pair<nda::array<double, 1>, nda::array<double, 2>> dos(tight_binding const &TB, int nkpts, int neps);

  std::pair<nda::array<double, 1>, nda::array<double, 1>> dos_patch(tight_binding const &TB, const nda::array<double, 2> &triangles, int neps,
                                                                    int ndiv);

 
} // namespace triqs::lattice