cyclat.hpp

Go to the documentation of this file.

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

 
#pragma once
 
#include "./mesh_iterator.hpp"
#include "./utils.hpp"
#include "../lattice/bravais_lattice.hpp"
#include "../utility/macros.hpp"
 
#include <fmt/ranges.h>
#include <h5/h5.hpp>
#include <nda/nda.hpp>
 
#include <array>
#include <cstdint>
#include <iostream>
#include <string>
 
namespace triqs::mesh {

  using lattice::bravais_lattice;


  class C2PY_RENAME(MeshCycLat) cyclat {
    public:
    using value_t = bravais_lattice::point_t;

    using index_t = std::array<long, 3>;

    using data_index_t = long;

    class C2PY_IGNORE mesh_point_t : public value_t {
      public:
      using mesh_t = cyclat;

      mesh_point_t() = default;

      mesh_point_t(std::array<long, 3> const &n, long d, uint64_t mhash, bravais_lattice const *bl_ptr)
         : value_t(n, bl_ptr), data_index_(d), mesh_hash_(mhash) {}

      [[nodiscard]] long data_index() const { return data_index_; }

      [[nodiscard]] value_t const &value() const { return *this; }

      [[nodiscard]] uint64_t mesh_hash() const noexcept { return mesh_hash_; }

      friend std::ostream &operator<<(std::ostream &sout, mesh_point_t const &mp) { return sout << mp.value(); }
 
      private:
      long data_index_    = 0;
      uint64_t mesh_hash_ = 0;
    };
 
    public:
    C2PY_DEPRECATED_PARAMETER_NAME(lattice : bl)
    cyclat(bravais_lattice const &bl, std::array<long, 3> const &dims)
       : bl_(bl),
         dims_(dims),
         size_(nda::stdutil::product(dims)),
         s2_(dims_[2]),
         s1_(dims_[1] * dims_[2]),
         units_(bl.units()),
         units_inv_(nda::linalg::inv(units_)),
         mesh_hash_(hash(nda::sum(bl.units()), dims[0], dims[1], dims[2])) {
      EXPECTS(dims_[0] > 0 and dims_[1] > 0 and dims_[2] > 0);
    }

    C2PY_IGNORE cyclat(bravais_lattice const &bl, nda::matrix<long> const &M) : cyclat(bl, std::array{M(0, 0), M(1, 1), M(2, 2)}) {
      EXPECTS((M.shape() == std::array{3l, 3l}));
      EXPECTS(nda::is_matrix_diagonal(M));
    }

    cyclat(bravais_lattice const &bl, long L) : cyclat{bl, std::array{L, (bl.ndim() >= 2 ? L : 1l), (bl.ndim() >= 3 ? L : 1)}} {}

    cyclat(long L1 = 1, long L2 = 1, long L3 = 1) : cyclat{bravais_lattice{nda::eye<double>(3)}, std::array{L1, L2, L3}} {}

    [[nodiscard]] bool is_index_valid(index_t const &n) const noexcept {
      for (auto i : nda::range(3))
        if (n[i] < 0 or n[i] >= dims_[i]) return false;
      return true;
    }

    [[nodiscard]] data_index_t to_data_index(index_t const &n) const {
      EXPECTS(is_index_valid(n));
      return n[0] * s1_ + n[1] * s2_ + n[2];
    }

    [[nodiscard]] C2PY_IGNORE data_index_t to_data_index(closest_mesh_point_t<value_t> const &cmp) const { return to_data_index(to_index(cmp)); }

    [[nodiscard]] index_t to_index(data_index_t d) const {
      EXPECTS(0 <= d and d < size());
      long const r0 = d % s1_;
      return {d / s1_, r0 / s2_, r0 % s2_};
    }

    [[nodiscard]] C2PY_IGNORE index_t to_index(closest_mesh_point_t<value_t> const &cmp) const { return cmp.value.index(); }

    [[nodiscard]] mesh_point_t operator[](long d) const { return {to_index(d), d, mesh_hash_, &bl_}; }

    [[nodiscard]] C2PY_IGNORE mesh_point_t operator[](closest_mesh_point_t<value_t> const &cmp) const { return (*this)[this->to_data_index(cmp)]; }

    [[nodiscard]] mesh_point_t operator()(index_t const &n) const { return {n, to_data_index(n), mesh_hash_, &bl_}; }

    [[nodiscard]] value_t to_value(index_t const &n) const {
      EXPECTS(is_index_valid(n));
      return {n, &bl_};
    }

    [[nodiscard]] C2PY_PROPERTY_GET(dims) auto const &dims() const { return dims_; }

    [[nodiscard]] C2PY_PROPERTY_GET(units) auto units() const { return nda::matrix_const_view<double>{units_}; }

    [[nodiscard]] C2PY_PROPERTY_GET(lattice) auto const &lattice() const noexcept { return bl_; }

    [[nodiscard]] C2PY_PROPERTY_GET(mesh_hash) uint64_t mesh_hash() const { return mesh_hash_; }

    [[nodiscard]] long size() const { return size_; }

    [[nodiscard]] index_t index_modulo(index_t const &n_tilde) const {
      return {positive_modulo(n_tilde[0], dims_[0]), positive_modulo(n_tilde[1], dims_[1]), positive_modulo(n_tilde[2], dims_[2])};
    }

    [[nodiscard]] auto begin() const { return mesh_iterator<cyclat>{.mesh_ptr = this, .data_index = 0}; }

    [[nodiscard]] auto cbegin() const { return begin(); }

    [[nodiscard]] auto end() const { return mesh_iterator<cyclat>{.mesh_ptr = this, .data_index = size()}; }

    [[nodiscard]] auto cend() const { return end(); }

    bool operator==(cyclat const &m) const { return bl_ == m.lattice() && dims_ == m.dims(); }

    friend std::ostream &operator<<(std::ostream &sout, cyclat const &m) {
      return sout << "Cyclic lattice mesh with linear dimensions " << m.dims() << " and an underlying " << m.lattice();
    }

    void serialize(auto &ar) const { ar & bl_ & dims_ & size_ & s2_ & s1_ & units_ & units_inv_ & mesh_hash_; }

    void deserialize(auto &ar) { ar & bl_ & dims_ & size_ & s2_ & s1_ & units_ & units_inv_ & mesh_hash_; }

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

    friend void h5_write(h5::group g, std::string const &name, cyclat const &m) {
      h5::group gr = g.create_group(name);
      h5::write_hdf5_format(gr, m); // NOLINT (downcasting to base class)
      h5::write(gr, "dims", m.dims_);
      h5::write(gr, "bravais_lattice", m.bl_);
    }

    friend void h5_read(h5::group g, std::string const &name, cyclat &m) {
      h5::group gr = g.open_group(name);
      h5::assert_hdf5_format(gr, m, true); // NOLINT (downcasting to base class)
 
      std::array<long, 3> dims{};
      if (gr.has_key("dims")) {
        h5::read(gr, "dims", dims);
      } else {
        // for backward compatibility
        auto M = h5::read<nda::matrix<long>>(gr, "periodization_matrix");
        dims   = {M(0, 0), M(1, 1), M(2, 2)};
      }
 
      bravais_lattice bl{};
      if (gr.has_key("bravais_lattice")) {
        h5::read(gr, "bravais_lattice", bl);
      } else {
        // for backward compatibility
        h5::read(gr, "bl", bl);
      }
 
      m = cyclat(bl, dims);
    }
 
    private:
    bravais_lattice bl_            = {};
    std::array<long, 3> dims_      = {0, 0, 0};
    long size_                     = 0;
    long s2_                       = 1;
    long s1_                       = 1;
    nda::matrix<double> units_     = nda::eye<double>(3);
    nda::matrix<double> units_inv_ = nda::eye<double>(3);
    uint64_t mesh_hash_            = 0;
  };

  auto evaluate(cyclat const &m, auto const &f, cyclat::index_t const &n_tilde) { return f(m.index_modulo(n_tilde)); }

  auto evaluate(cyclat const &m, auto const &f, cyclat::value_t const &r_n_tilde) { return evaluate(m, f, r_n_tilde.index()); }

 
} // namespace triqs::mesh