imperative_operator.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-2023 Simons Foundation
// Copyright (c) 2016 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: Maxime Charlebois, Michel Ferrero, Igor Krivenko, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "./fundamental_operator_set.hpp"
#include "./hilbert_space.hpp"
#include "../operators/many_body_operator.hpp"
 
#include <algorithm>
#include <cstdint>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace triqs::hilbert_space {


  template <typename HS, typename T = double, bool UseMap = false> class imperative_operator {
    private:
    // Precomputed bitmask representation of a single monomial.
    struct one_term_t {
      T coeff;
      uint64_t d_mask, dag_mask, d_count_mask, dag_count_mask;
    };
 
    public:
    using coeff_t = T;

    using hilbert_map_t = std::vector<int>;

    imperative_operator() = default;

    imperative_operator(triqs::operators::many_body_operator_generic<coeff_t> const &op, fundamental_operator_set const &fops,
                        hilbert_map_t hmap = hilbert_map_t(), std::vector<sub_hilbert_space> const *subspaces = nullptr)
       : subspaces_(subspaces), hmap_(std::move(hmap)) {
 
      if ((hmap_.size() == 0) != !UseMap) TRIQS_RUNTIME_ERROR << "Internal error";
 
      // ordering matching the canonical form
      auto greater = [&fops](triqs::operators::canonical_ops_t const &op1, triqs::operators::canonical_ops_t const &op2) {
        if (op1.dagger != op2.dagger) return op2.dagger;
        return op1.dagger ? (fops[op1.indices] > fops[op2.indices]) : (fops[op1.indices] < fops[op2.indices]);
      };
 
      // The goal here is to have a transcription of the many_body_operator in terms
      // of simple vectors (maybe the code below could be more elegant)
      for (auto const &term : op) {
        // canonical-order working copy; each swap flips the fermionic sign of the coefficient
        auto monomial = term.monomial;
        auto coef     = term.coef;
 
        // Sort monomial according to the order established by fops
        int n        = monomial.size();
        bool swapped = true;
        while (swapped) {
          swapped = false;
          for (int i = 1; i < n; ++i) {
            if (greater(monomial[i - 1], monomial[i])) {
              using std::swap;
              swap(monomial[i - 1], monomial[i]);
              swapped = true;
              coef *= coeff_t(-1);
            }
          }
          --n;
        }
 
        // Given the environment variable CHECK_ISSUE819 was set by the user
        // throw an exception if the result of this model was effected by issue 819
        // https://github.com/TRIQS/triqs/issues/819
        static const bool check_issue819 = std::getenv("CHECK_ISSUE819");
        if (check_issue819 && term.coef != coef)
          TRIQS_RUNTIME_ERROR << "ERROR: The Atom-Diag result of this model is affected by issue 819 (https://github.com/TRIQS/triqs/issues/819).\n"
                                 "If you have solved the same model with release 2.2.0, 2.2.1 or 3.0.0 of TRIQS the result was incorrect.";
 
        // build bitmask representation of the canonical monomial
        std::vector<int> dag, ndag;
        uint64_t d_mask = 0, dag_mask = 0;
        for (auto const &canonical_op : monomial) {
          (canonical_op.dagger ? dag : ndag).push_back(fops[canonical_op.indices]);
          (canonical_op.dagger ? dag_mask : d_mask) |= (uint64_t(1) << fops[canonical_op.indices]);
        }
        // parity mask: bit i set iff i is not in d AND the number of d-bits > i is odd
        // used to read off the fermionic sign via parity_number_of_bits in operator()
        auto compute_count_mask = [](std::vector<int> const &d) {
          uint64_t mask = 0;
          bool is_on    = (d.size() % 2 == 1);
          for (int i = 0; i < 64; ++i) {
            if (std::ranges::find(d, i) != d.end())
              is_on = !is_on;
            else if (is_on)
              mask |= (uint64_t(1) << i);
          }
          return mask;
        };
        uint64_t d_count_mask = compute_count_mask(ndag), dag_count_mask = compute_count_mask(dag);
        terms_.push_back(one_term_t{coeff_t(coef), d_mask, dag_mask, d_count_mask, dag_count_mask});
      }
    }

    template <typename F> void update(F f) {
      for (auto &t_i : terms_) f(t_i.coeff);
    }

    [[nodiscard]] bool is_empty() const { return (terms_.size() == 0); }

    template <typename S, typename... Args> S operator()(S const &psi, Args &&...args) const {
 
      S target_st    = get_target_st(psi);
      auto const &hs = psi.get_hilbert();
 
      using amplitude_t = typename S::value_type;
 
      for (int i = 0; i < terms_.size(); ++i) { // loop over monomials
        auto M = terms_[i];
        foreach (psi, [M, &target_st, hs, ... args = std::forward<Args>(args)](int j, typename S::value_type amplitude) {
          fock_state_t f2 = hs.get_fock_state(j);
          if ((f2 & M.d_mask) != M.d_mask) return;
          f2 &= ~M.d_mask;
          if (((f2 ^ M.dag_mask) & M.dag_mask) != M.dag_mask) return;
          fock_state_t f3    = ~(~f2 & ~M.dag_mask);
          auto sign_is_minus = parity_number_of_bits((f2 & M.d_count_mask) ^ (f3 & M.dag_count_mask));
          // update state vector in target Hilbert space
          auto ind = target_st.get_hilbert().get_state_index(f3);
          target_st(ind) += amplitude * apply_if_possible(M.coeff, args...) * (sign_is_minus ? -amplitude_t(1) : amplitude_t(1));
        }); // foreach
      }
      return target_st;
    }
 
    private:
    // Return a zero state in the target Hilbert space of |psi> under this operator, or a default-constructed (empty)
    // state when UseMap=true and the connection map sends |psi>'s subspace to -1.
    template <typename S> S get_target_st(S const &psi) const {
      if constexpr (UseMap) {
        auto n = hmap_[psi.get_hilbert().get_index()];
        if (n == -1) return S{};
        return S{(*subspaces_)[n]};
      } else {
        return S(psi.get_hilbert());
      }
    }
 
    // Return true if the number of set bits in v is odd, false if it is even.
    static bool parity_number_of_bits(uint64_t v) {
      // http://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetNaive
      // v ^= v >> 16;
      // only ok until 16 orbitals ! assert this or put the >> 16
      v ^= v >> 8;
      v ^= v >> 4;
      v ^= v >> 2;
      v ^= v >> 1;
      return v & 0x01;
    }
 
    // Return x(args...) when args... is non-empty (coeff_t must be callable with them), x otherwise.
    // Supports parameterized operators whose coefficients are themselves callable objects.
    template <typename... Args> static auto apply_if_possible(coeff_t const &x, Args &&...args) -> std::invoke_result_t<coeff_t, Args...> {
      return x(std::forward<Args>(args)...);
    }
    static auto apply_if_possible(coeff_t const &x) -> coeff_t { return x; }
 
    private:
    std::vector<one_term_t> terms_;
    std::vector<sub_hilbert_space> const *subspaces_ = nullptr;
    hilbert_map_t hmap_;
  };

 
} // namespace triqs::hilbert_space