extractors.hpp

Go to the documentation of this file.

// Copyright (c) 2015-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2015-2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2021 Simons Foundation
// Copyright (c) 2015-2017 Igor Krivenko
//
// 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: Thomas Ayral, Michel Ferrero, Igor Krivenko, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "../many_body_operator.hpp"
#include "../../arrays.hpp"
#include "../../hilbert_space/fundamental_operator_set.hpp"
#include "../../utility/first_include.hpp"
 
#include <complex>
#include <map>
#include <string>
#include <tuple>
#include <variant>
 
namespace triqs::operators::utils {

 
  // Elevate `nda::array` to the `triqs::operators::utils` namespace.
  using nda::array;
 
  // Elevate triqs::hilbert_space::fundamental_operator_set::indices_t to the `triqs::operators::utils` namespace.
  using indices_t = hilbert_space::fundamental_operator_set::indices_t;

  template <typename T> using op_t = operators::many_body_operator_generic<T>;

  template <typename T> using dict2_t = std::map<std::tuple<indices_t, indices_t>, T>;

  template <typename T> using dict4_t = std::map<std::tuple<indices_t, indices_t, indices_t, indices_t>, T>;

  template <int N> using real_or_complex_array = std::variant<array<double, N>, array<std::complex<double>, N>>;

  template <typename T> using block_matrix_t = array<nda::matrix<T>, 1>;

  template <typename T> [[nodiscard]] dict2_t<T> extract_h_dict(op_t<T> const &h, bool ignore_irrelevant = false) {
    auto h_dict = dict2_t<T>{};
 
    for (auto const &term : h) {
      auto const &coef = term.coef;
      auto const &m    = term.monomial;
 
      if (m.size() == 2) {
        if (!(m[0].dagger && !m[1].dagger)) {
          if (!ignore_irrelevant) TRIQS_RUNTIME_ERROR << "extract_h_dict: monomial is not of the form C^+(i) C(j)";
        } else { // everything ok
          h_dict.insert({std::make_tuple(m[0].indices, m[1].indices), coef});
        }
      } else {
        if (!ignore_irrelevant) TRIQS_RUNTIME_ERROR << "extract_h_dict: monomial must have 2 operators";
      }
    }
 
    return h_dict;
  }

  template <typename T> [[nodiscard]] dict2_t<T> extract_U_dict2(op_t<T> const &h, bool ignore_irrelevant = false) {
    auto U_dict = dict2_t<T>{};
 
    for (auto const &term : h) {
      auto const &coef = term.coef;
      auto const &m    = term.monomial;
 
      if (m.size() == 4) {
        if (!(m[0].dagger && m[1].dagger && !m[2].dagger && !m[3].dagger) || (m[0].indices != m[3].indices) || (m[1].indices != m[2].indices)) {
          if (!ignore_irrelevant) TRIQS_RUNTIME_ERROR << "extract_U_dict2: monomial is not of the form C^+(i) C^+(j) C(j) C(i)";
        } else { //everything ok
          U_dict.insert({std::make_tuple(m[0].indices, m[1].indices), coef});
          U_dict.insert({std::make_tuple(m[1].indices, m[0].indices), coef});
        }
      } else {
        if (!ignore_irrelevant) TRIQS_RUNTIME_ERROR << "extract_U_dict2: monomial must have 4 operators";
      }
    }
 
    return U_dict;
  }

  template <typename T> [[nodiscard]] dict4_t<T> extract_U_dict4(op_t<T> const &h, bool ignore_irrelevant = false) {
    auto U_dict = dict4_t<T>{};
 
    for (auto const &term : h) {
      T const &coef = term.coef;
      auto const &m = term.monomial;
 
      if (m.size() == 4) {
        if (!(m[0].dagger && m[1].dagger && !m[2].dagger && !m[3].dagger)) {
          if (!ignore_irrelevant) TRIQS_RUNTIME_ERROR << "extract_U_dict4: monomial is not of the form C^+(i) C^+(j) C(l) C(k)";
        } else { // everything ok
          U_dict.insert({std::make_tuple(m[0].indices, m[1].indices, m[3].indices, m[2].indices), T(0.5) * coef});
          U_dict.insert({std::make_tuple(m[1].indices, m[0].indices, m[2].indices, m[3].indices), T(0.5) * coef});
          U_dict.insert({std::make_tuple(m[0].indices, m[1].indices, m[2].indices, m[3].indices), T(-0.5) * coef});
          U_dict.insert({std::make_tuple(m[1].indices, m[0].indices, m[3].indices, m[2].indices), T(-0.5) * coef});
        }
      } else {
        if (!ignore_irrelevant) TRIQS_RUNTIME_ERROR << "extract_U_dict4: monomial must have 4 operators";
      }
    }
 
    return U_dict;
  }

  template <typename T = double, typename D, std::size_t R = std::tuple_size_v<typename D::key_type>>
  [[nodiscard]] array<T, R> dict_to_matrix(D const &dict, hilbert_space::fundamental_operator_set const &fs) {
    using triqs::hilbert_space::format_indices;
 
    // check if a given alpha_i is in the fops and return its linear index
    auto indices_to_linear = [&fs](indices_t const &indices) {
      if (!fs.has_indices(indices))
        TRIQS_RUNTIME_ERROR << "dict_to_matrix: key [" << format_indices(indices) << "] of dict not in fundamental_operator_set/gf_struct";
      return fs[indices];
    };
 
    // create and fill the coefficient array
    auto arr = nda::zeros<T>(nda::stdutil::make_initialized_array<R>(fs.size()));
    for (auto const &[key, value] : dict) std::apply([&](auto const &...ks) -> auto & { return arr(indices_to_linear(ks)...); }, key) = T(value);
    return arr;
  }

  template <typename D, std::size_t R = std::tuple_size_v<typename D::key_type>>
  [[nodiscard]] real_or_complex_array<R> dict_to_variant_matrix(D const &dict, hilbert_space::fundamental_operator_set const &fs) {
    for (auto const &[key, value] : dict) {
      if (!value.is_real()) return dict_to_matrix<std::complex<double>>(dict, fs);
    }
    return dict_to_matrix<double>(dict, fs);
  }

  template <typename T> [[nodiscard]] op_t<T> filter_op(op_t<T> const &h, long len) {
    auto h_filtered = op_t<T>{};
    for (auto const &term : h) {
      if (term.monomial.size() == len) h_filtered += term;
    }
    return h_filtered;
  }

  template <typename T> [[nodiscard]] op_t<T> quadratic_terms(op_t<T> const &h) { return filter_op(h, 2); }

  template <typename T> [[nodiscard]] op_t<T> quartic_terms(op_t<T> const &h) { return filter_op(h, 4); }

  template <typename T>
  [[nodiscard]] block_matrix_t<T> block_matrix_from_op(op_t<T> const &h, hilbert_space::gf_struct_t const &gf_struct,
                                                       bool ignore_irrelevant = false) {
    auto n_bl   = gf_struct.size();
    auto bl_mat = nda::array<nda::matrix<T>, 1>(n_bl);
    for (auto bl : range(n_bl)) {
      auto bl_size = gf_struct[bl].second;
      bl_mat[bl]   = nda::zeros<T>(bl_size, bl_size);
    }
 
    auto name_to_bl = [&](auto &bl_name) {
      auto it =
         std::find_if(cbegin(gf_struct), cend(gf_struct), [&bl_name](auto const &blname_and_size) { return bl_name == blname_and_size.first; });
      return std::distance(cbegin(gf_struct), it);
    };
 
    for (auto const &term : h) {
      auto const &m = term.monomial;
      T const &coef = term.coef;
 
      if (m.size() == 2 and m[0].dagger and not m[1].dagger and m[0].indices[0] == m[1].indices[0] and m[0].indices.size() == 2
          and m[1].indices.size() == 2) {
        auto bl_name = std::get<std::string>(m[0].indices[0]);
        auto op1_idx = std::get<long>(m[0].indices[1]);
        auto op2_idx = std::get<long>(m[1].indices[1]);
 
        bl_mat[name_to_bl(bl_name)](op1_idx, op2_idx) = coef;
      } else {
        if (!ignore_irrelevant) TRIQS_RUNTIME_ERROR << "block_matrix_from_op: Operator term is not of the form 'coeff * c_dag{bl,i} * c_{bl,j}'";
      }
    }
 
    return bl_mat;
  }

  template <typename T> [[nodiscard]] op_t<T> op_from_block_matrix(block_matrix_t<T> const &bl_mat, hilbert_space::gf_struct_t const &gf_struct) {
    EXPECTS(bl_mat.size() == gf_struct.size());
    auto h = op_t<T>{};
    for (long bl = 0; auto const &[bl_name, bl_size] : gf_struct) {
      EXPECTS(bl_mat[bl].shape() == (std::array{bl_size, bl_size}));
      for (auto [i, j] : itertools::product_range(bl_size, bl_size)) h += bl_mat[bl](i, j) * c_dag(bl_name, i) * c(bl_name, j);
      ++bl;
    }
    return h;
  }

 
} // namespace triqs::operators::utils