state.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) 2015-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: Michel Ferrero, Igor Krivenko, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "./hilbert_space.hpp"
#include "../arrays.hpp"
#include "../utility/numeric_ops.hpp"
 
#include <boost/operators.hpp>
 
#include <cstddef>
#include <iostream>
#include <unordered_map>
 
namespace triqs::hilbert_space {

 
  // Forward declaration.
  template <typename HS, typename T, bool BasedOnMap> class state {};

  template <typename HS, typename T, bool BasedOnMap> auto make_zero_state(state<HS, T, BasedOnMap> const &phi) {
    return state<HS, T, BasedOnMap>{phi.get_hilbert()};
  }

  template <typename HS, typename T> class state<HS, T, true> : boost::additive<state<HS, T, true>>, boost::multiplicative<state<HS, T, true>, T> {
    public:
    using value_type = T;

    using hilbert_space_t = HS;

    using amplitude_t = std::unordered_map<fock_state_t, value_type>;

    state() : hs_ptr_(nullptr) {}

    state(hilbert_space_t const &hs) : hs_ptr_(&hs) {}

    state(hilbert_space_t const &hs, fock_state_t f) : hs_ptr_(&hs) { map_[f] = value_type(1.0); }

    [[nodiscard]] int size() const { return hs_ptr_->size(); }

    [[nodiscard]] int nterms() const { return map_.size(); }

    value_type &operator()(fock_state_t f) { return map_[f]; }

    value_type const &operator()(fock_state_t f) const { return map_[f]; }

    state &operator+=(state const &phi) {
      for (auto const &[f, b_f] : phi.map_) {
        auto [it, inserted] = map_.insert({f, b_f});
        if (!inserted) it->second += b_f;
      }
      prune();
      return *this;
    }

    state &operator-=(state const &phi) {
      for (auto const &[f, b_f] : phi.map_) {
        auto [it, inserted] = map_.insert({f, -b_f});
        if (!inserted) it->second -= b_f;
      }
      prune();
      return *this;
    }

    state &operator*=(value_type x) {
      for (auto &[f, a_f] : map_) { a_f *= x; }
      prune();
      return *this;
    }

    state &operator/=(value_type x) { return operator*=(value_type(1.0) / x); }

    friend value_type dot_product(state const &phi, state const &psi) {
      value_type res(0.0);
      for (auto const &[f, b_f] : phi.map_) {
        using triqs::utility::conj;
        if (psi.map_.count(f) == 1) res += conj(b_f) * psi.map_.at(f);
      }
      return res;
    }

    template <typename F> friend void foreach (state const &psi, F lambda) {
      const_cast<state &>(psi).prune(); // NOLINT (why do we const cast here?)
      for (auto const &[f, a_f] : psi.map_) lambda(f, a_f);
    }

    hilbert_space_t const &get_hilbert() const { return *hs_ptr_; }

    void set_hilbert(hilbert_space_t const &new_hs) { hs_ptr_ = &new_hs; }
 
    private:
    // Discard all amplitudes smaller than a specific threshold.
    void prune() {
      for (auto it = map_.begin(); it != map_.end();) {
        using triqs::utility::is_zero;
        if (is_zero(it->second)) {
          it = map_.erase(it);
        } else {
          ++it;
        }
      }
    }
 
    private:
    const hilbert_space_t *hs_ptr_;
    amplitude_t map_;
  };

  template <typename HS, typename T> class state<HS, T, false> : boost::additive<state<HS, T, false>>, boost::multiplicative<state<HS, T, false>, T> {
    public:
    using value_type = T;

    using hilbert_space_t = HS;

    using amplitude_t = nda::vector<value_type>;

    state() : hs_ptr_(nullptr) {}

    state(hilbert_space_t const &hs) : hs_ptr_(&hs), vec_(amplitude_t::zeros(hs.size())) {}

    state(hilbert_space_t const &hs, fock_state_t f) : hs_ptr_(&hs), vec_(amplitude_t::zeros(hs.size())) { vec_[f] = value_type(1.0); }

    [[nodiscard]] int size() const { return hs_ptr_->size(); }

    value_type &operator()(int i) { return vec_[i]; }

    value_type const &operator()(int i) const { return vec_[i]; }

    state &operator+=(state const &phi) {
      vec_ += phi.vec_;
      return *this;
    }

    state &operator-=(state const &phi) {
      vec_ -= phi.vec_;
      return *this;
    }

    state &operator*=(value_type x) {
      vec_ *= x;
      return *this;
    }

    state &operator/=(value_type x) {
      vec_ /= x;
      return *this;
    }

    friend value_type dot_product(state const &phi, state const &psi) { return nda::blas::dotc(phi.vec_, psi.vec_); }

    template <typename F> friend void foreach (state const &psi, F lambda) {
      const auto dim = psi.size();
      for (size_t i = 0; i < dim; ++i) lambda(i, psi(i));
    }

    amplitude_t const &amplitudes() const { return vec_; }

    amplitude_t &amplitudes() { return vec_; }

    hilbert_space_t const &get_hilbert() const { return *hs_ptr_; }

    void set_hilbert(hilbert_space_t const &new_hs) { hs_ptr_ = &new_hs; }
 
    private:
    const hilbert_space_t *hs_ptr_;
    amplitude_t vec_;
  };

  template <typename HS, typename T, bool BasedOnMap> std::ostream &operator<<(std::ostream &sout, state<HS, T, BasedOnMap> const &psi) {
    bool something_written = false;
    auto const &hs         = psi.get_hilbert();
 
    foreach (psi, [&sout, hs, &something_written](int i, auto a_f_i) {
      using triqs::utility::is_zero;
      if (!is_zero(a_f_i)) {
        sout << " +(" << a_f_i << ")"
             << "|" << hs.get_fock_state(i) << ">";
        something_written = true;
      }
    });
 
    if (!something_written) sout << 0;
    return sout;
  }

  template <typename TargetState, typename OriginalState> TargetState project(OriginalState const &psi, hilbert_space const &proj_hs) {
    TargetState proj_psi(proj_hs);
    auto const &hs = psi.get_hilbert();
    foreach (psi, [&](int i, auto a_f_i) { proj_psi(hs.get_fock_state(i)) = a_f_i; });
    return proj_psi;
  }

  template <typename TargetState, typename OriginalState> TargetState project(OriginalState const &psi, sub_hilbert_space const &proj_hs) {
    TargetState proj_psi(proj_hs);
    auto const &hs = psi.get_hilbert();
    foreach (psi, [&](int i, auto a_f_i) {
      auto f = hs.get_fock_state(i);
      if (proj_hs.has_state(f)) proj_psi(proj_hs.get_state_index(f)) = a_f_i;
    });
    return proj_psi;
  }

 
} // namespace triqs::hilbert_space