atom_diag.hpp

Go to the documentation of this file.

// Copyright (c) 2017-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2017-2018 Centre national de la recherche scientifique (CNRS)
// Copyright (c) 2018-2022 Simons Foundation
// Copyright (c) 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: Maxime Charlebois, Michel Ferrero, Igor Krivenko, Olivier Parcollet, Hugo U. R. Strand, Nils Wentzell

 
#pragma once
 
#include "../arrays.hpp"
#include "../hilbert_space/fundamental_operator_set.hpp"
#include "../hilbert_space/hilbert_space.hpp"
#include "../operators/many_body_operator.hpp"
#include "../utility/macros.hpp"
 
#include <h5/h5.hpp>
#include <itertools/itertools.hpp>
#include <nda/h5.hpp>
#include <nda/nda.hpp>
 
#include <complex>
#include <initializer_list>
#include <iostream>
#include <string>
#include <type_traits>
#include <vector>
 
namespace triqs::atom_diag {
 
  // Import tools from the hilbert_space and arrays namespace.
  using namespace triqs::hilbert_space;
  using namespace triqs::arrays;
 
  // Forward declarations.
  template <bool C> class atom_diag;
  template <bool C> struct atom_diag_worker;
  template <bool C> std::ostream &operator<<(std::ostream &, atom_diag<C> const &);
  template <bool C> void h5_write(h5::group, std::string const &, atom_diag<C> const &);
  template <bool C> void h5_read(h5::group, std::string const &, atom_diag<C> &);


  using indices_t = fundamental_operator_set::indices_t;

  using quantum_number_t = double;

  template <bool Complex> class atom_diag {
    public:
    using scalar_t = std::conditional_t<Complex, std::complex<double>, double>;

    using matrix_t = matrix<scalar_t>;

    using block_matrix_t = std::vector<matrix_t>;

    using full_hilbert_space_state_t = vector<scalar_t>;

    using many_body_op_t = triqs::operators::many_body_operator_generic<scalar_t>;

    using fops_data_t = triqs::hilbert_space::fundamental_operator_set::data_t;

    struct C2PY_IGNORE eigensystem_t {
      vector<double> eigenvalues;

      matrix_t unitary_matrix;

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

      C2PY_IGNORE friend void mpi_broadcast(eigensystem_t &es, mpi::communicator c = {}, int root = 0) {
        mpi::broadcast(es.eigenvalues, c, root);
        mpi::broadcast(es.unitary_matrix, c, root);
      }

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

      friend void h5_write(h5::group g, std::string const &name, eigensystem_t const &es) {
        auto gr = g.create_group(name);
        h5_write(gr, "eigenvalues", es.eigenvalues);
        h5_write(gr, "unitary_matrix", es.unitary_matrix);
      }

      friend void h5_read(h5::group g, std::string const &name, eigensystem_t &es) {
        auto gr = g.open_group(name);
        h5_read(gr, "eigenvalues", es.eigenvalues);
        h5_read(gr, "unitary_matrix", es.unitary_matrix);
      }
    };

    struct C2PY_IGNORE op_block_mat_t {
      array<long, 1> connection;

      std::vector<matrix_t> block_mat;

      op_block_mat_t(int n_blocks) : connection(n_blocks), block_mat(n_blocks) { connection(range::all) = -1; };

      int n_blocks() const { return block_mat.size(); }

      friend std::ostream &operator<<(std::ostream &os, op_block_mat_t const &op_mat) {
        os << "Operator block matrix:\n"
           << " n_blocks = " << op_mat.n_blocks() << "\n";
        for (auto bidx : range(op_mat.n_blocks())) {
          auto &m = op_mat.block_mat[bidx];
          auto b2 = op_mat.connection(bidx);
          os << "Block: "
             << "blocks (" << bidx << ", " << b2 << ") "
             << "size (" << first_dim(m) << ", " << second_dim(m) << ") " << op_mat.block_mat[bidx] << "\n";
        }
        return os;
      };
    };

    C2PY_IGNORE atom_diag() = default;

    atom_diag(many_body_op_t const &h, fundamental_operator_set const &fops);

    atom_diag(many_body_op_t const &h, fundamental_operator_set const &fops, many_body_op_t const &hyb);

    atom_diag(many_body_op_t const &h, fundamental_operator_set const &fops, int n_min, int n_max);

    atom_diag(many_body_op_t const &h, fundamental_operator_set const &fops, std::vector<many_body_op_t> const &qn_vector);

    C2PY_IGNORE atom_diag(many_body_op_t const &h, fundamental_operator_set const &fops, std::initializer_list<many_body_op_t> const &init_lst)
       : atom_diag(h, fops, std::vector<many_body_op_t>{init_lst}) {};

    C2PY_IGNORE many_body_op_t const &get_h_atomic() const { return h_atomic; }

    C2PY_PROPERTY_GET(h_atomic) operators::many_body_operator get_h_atomic_as_mbop() const { return h_atomic; }

    C2PY_IGNORE fundamental_operator_set const &get_fops() const { return fops; }

    C2PY_PROPERTY_GET(fops) fops_data_t const &get_fops_as_data() const { return fops.data(); }

    C2PY_IGNORE class hilbert_space const &get_full_hilbert_space() const { return full_hs; }

    C2PY_PROPERTY_GET(full_hilbert_space_dim) int get_full_hilbert_space_dim() const { return full_hs.size(); }

    C2PY_PROPERTY_GET(n_subspaces) int n_subspaces() const { return eigensystems.size(); }

    int get_subspace_dim(int sp_index) const { return eigensystems[sp_index].eigenvalues.size(); }

    std::vector<int> get_subspace_dims() const {
      auto dims = std::vector<int>(n_subspaces());
      for (long i : range(n_subspaces())) dims[i] = get_subspace_dim(i);
      return dims;
    }

    C2PY_IGNORE std::vector<fock_state_t> const &get_fock_states(int sp_index) const { return sub_hilbert_spaces[sp_index].get_all_fock_states(); }

    C2PY_PROPERTY_GET(fock_states) std::vector<std::vector<fock_state_t>> get_fock_states() const {
      std::vector<std::vector<fock_state_t>> fock_states(n_subspaces());
      for (auto i : range(n_subspaces())) fock_states[i] = sub_hilbert_spaces[i].get_all_fock_states();
      return fock_states;
    }

    matrix<scalar_t> const &get_unitary_matrix(int sp_index) const { return eigensystems[sp_index].unitary_matrix; }

    C2PY_PROPERTY_GET(unitary_matrices) std::vector<matrix<scalar_t>> get_unitary_matrices() const {
      std::vector<matrix<scalar_t>> umat(n_subspaces());
      for (auto i : range(n_subspaces())) umat[i] = get_eigensystems()[i].unitary_matrix;
      return umat;
    }

    int flatten_subspace_index(int sp_index, int i) const { return first_eigenstate_of_subspace[sp_index] + i; }

    C2PY_IGNORE range index_range_of_subspace(int sp_index) const {
      return range{first_eigenstate_of_subspace[sp_index], first_eigenstate_of_subspace[sp_index] + get_subspace_dim(sp_index)};
    }

    C2PY_IGNORE std::vector<eigensystem_t> const &get_eigensystems() const { return eigensystems; }

    double get_eigenvalue(int sp_index, int i) const { return eigensystems[sp_index].eigenvalues[i]; }

    C2PY_PROPERTY_GET(energies) std::vector<std::vector<double>> get_energies() const;

    C2PY_PROPERTY_GET(quantum_numbers) std::vector<std::vector<quantum_number_t>> const &get_quantum_numbers() const { return quantum_numbers; }

    C2PY_PROPERTY_GET(gs_energy) double get_gs_energy() const { return gs_energy; }

    C2PY_PROPERTY_GET(vacuum_subspace_index) long get_vacuum_subspace_index() const { return vacuum_subspace_index; }

    C2PY_PROPERTY_GET(vacuum_state) full_hilbert_space_state_t const &get_vacuum_state() const { return vacuum; }

    long c_connection(int op_linear_index, int sp_index) const { return annihilation_connection(op_linear_index, sp_index); }

    long cdag_connection(int op_linear_index, int sp_index) const { return creation_connection(op_linear_index, sp_index); }

    matrix_t const &c_matrix(int op_linear_index, int sp_index) const { return c_matrices[op_linear_index][sp_index]; }

    matrix_t const &cdag_matrix(int op_linear_index, int sp_index) const { return cdag_matrices[op_linear_index][sp_index]; }

    C2PY_IGNORE std::pair<int, matrix_t> get_matrix_element_of_monomial(operators::monomial_t const &op_vec, int B) const;

    C2PY_IGNORE op_block_mat_t get_op_mat(many_body_op_t const &op) const;

    bool operator==(atom_diag const &rhs) const = default;

    [[nodiscard]] static std::string hdf5_format() { return Complex ? "AtomDiagComplex" : "AtomDiagReal"; }
 
    // Friend declarations.
    friend struct atom_diag_worker<Complex>;
    friend std::ostream &operator<< <Complex>(std::ostream &os, atom_diag const &ss);
    friend void h5_write<Complex>(h5::group gr, std::string const &name, atom_diag const &);
    friend void h5_read<Complex>(h5::group gr, std::string const &name, atom_diag &);
 
    private:
    void fill_first_eigenstate_of_subspace();
    void compute_vacuum();
 
    private:
    many_body_op_t h_atomic;                                    // Hamiltonian
    fundamental_operator_set fops;                              // Keep it to compute the Green's function
    class hilbert_space full_hs;                                // Full Hilbert space of the problem
    std::vector<sub_hilbert_space> sub_hilbert_spaces;          // The invariant subspaces, i.e. the lists of Fock states
    std::vector<eigensystem_t> eigensystems;                    // Eigensystem in each subspace
    matrix<long> creation_connection;                           // creation_connection(i, B) -> B', target subspace of c†_i acting on B
    matrix<long> annihilation_connection;                       // idem for annihilation operators
    std::vector<std::vector<matrix_t>> cdag_matrices;           // cdag_matrices[i][B] = block of c†_i mapping B -> B'
    std::vector<std::vector<matrix_t>> c_matrices;              // idem for annihilation operators
    double gs_energy{0.0};                                      // Energy of the ground state
    long vacuum_subspace_index{-1};                             // Invariant subspace containing |0>
    full_hilbert_space_state_t vacuum;                          // Vacuum vector (in the eigenbasis)
    std::vector<std::vector<quantum_number_t>> quantum_numbers; // Values of the quantum numbers for each subspace
    std::vector<int> first_eigenstate_of_subspace;              // Index of the first eigenstate of each subspace
  };

 
} // namespace triqs::atom_diag