bravais_lattice.hpp

Go to the documentation of this file.

// Copyright (c) 2014-2018 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
// Copyright (c) 2014-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: Michel Ferrero, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "../arrays.hpp"
#include "../utility/macros.hpp"
 
#include <fmt/ranges.h>
#include <h5/h5.hpp>
#include <nda/nda.hpp>
 
#include <array>
#include <iostream>
#include <optional>
#include <ranges>
#include <string>
#include <vector>
 
namespace triqs::lattice {


  using r_t = nda::vector<double>;

  using matrix_t = nda::matrix<double>;

  class bravais_lattice {
    public:
    using index_t = std::array<long, 3>;

    class C2PY_IGNORE point_t {
      public:
      point_t() = default;

      point_t(index_t const &n, bravais_lattice const *bl_ptr) : index_(n), bl_ptr_(bl_ptr) {}

      [[nodiscard]] auto index() const { return index_; }

      [[nodiscard]] auto const &lattice() const { return *bl_ptr_; }

      explicit operator r_t() const {
        if (rval_)
          return *rval_;
        else
          return *(rval_ = bl_ptr_->lattice_to_real_coordinates(index_));
      }

      [[deprecated("() is deprecated for a cyclat::mesh_point_t. Use [] instead")]] double operator()(int d) const { return r_t(*this)[d]; }

      double operator[](int i) const { return r_t(*this)[i]; }

      [[nodiscard]] point_t operator+(point_t const &pt) const {
        EXPECTS(*bl_ptr_ == *(pt.bl_ptr_));
        return {index_ + pt.index_, bl_ptr_};
      }

      [[nodiscard]] point_t operator-(point_t const &pt) const {
        EXPECTS(*bl_ptr_ == *(pt.bl_ptr_));
        return {index_ - pt.index_, bl_ptr_};
      }

      [[nodiscard]] point_t operator-() const { return {-index_, bl_ptr_}; }

      friend std::ostream &operator<<(std::ostream &sout, point_t const &pt) { return sout << static_cast<r_t>(pt); }
 
      private:
      index_t index_                   = {0, 0, 0};
      bravais_lattice const *bl_ptr_   = nullptr;
      mutable std::optional<r_t> rval_ = {};
    };

    bravais_lattice() : bravais_lattice(nda::eye<double>(3)) {}

    bravais_lattice(matrix_t const &units, std::vector<r_t> orbital_positions = std::vector<r_t>{{0, 0, 0}},
                    std::vector<std::string> atom_orb_name = {});

    [[nodiscard]] bool contains(r_t const &) const { return true; }

    [[nodiscard]] int ndim() const { return ndim_; }

    [[nodiscard]] matrix_t const &units() const { return units_; }

    [[nodiscard]] long n_orbitals() const { return static_cast<long>(atom_orb_pos_.size()); }

    [[nodiscard]] auto const &orbital_positions() const { return atom_orb_pos_; }

    [[nodiscard]] auto const &orbital_names() const { return atom_orb_name_; }

    template <typename R> [[nodiscard]] r_t lattice_to_real_coordinates(R const &v) const {
      return nda::transpose(units_)(nda::range::all, nda::range(ndim())) * nda::basic_array_view{v}(nda::range(ndim()));
    }

    template <typename R> [[nodiscard]] r_t real_to_lattice_coordinates(R const &v_tilde) const {
      return nda::transpose(units_inv_)(nda::range::all, nda::range(ndim())) * nda::basic_array_view{v_tilde}(nda::range(ndim()));
    }

    [[nodiscard]] C2PY_IGNORE point_t to_point(index_t const &n) const { return {n, this}; }

    bool operator==(bravais_lattice const &bl) const {
      return units_ == bl.units() && ndim_ == bl.ndim() && atom_orb_pos_ == bl.orbital_positions() && atom_orb_name_ == bl.orbital_names();
    }

    friend std::ostream &operator<<(std::ostream &sout, bravais_lattice const &bl) {
      auto str = fmt::format("Bravais Lattice in {} dimensions with {} orbital(s):\n", bl.ndim(), bl.n_orbitals());
      str += fmt::format("  Basis vectors: {}\n",
                         std::views::transform(nda::range(bl.ndim()), [&bl](int i) { return bl.units()(i, nda::range(bl.ndim())); }));
      str += fmt::format("  Orbital positions: {}\n", bl.orbital_positions());
      str += fmt::format("  Orbital names: {}\n", bl.orbital_names());
      return sout << str;
    }

    void serialize(auto &ar) const { ar & ndim_ & units_ & atom_orb_pos_ & atom_orb_name_; }

    void deserialize(auto &ar) { ar & ndim_ & units_ & atom_orb_pos_ & atom_orb_name_; }

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

    friend void h5_write(h5::group g, std::string const &name, bravais_lattice const &bl);

    friend void h5_read(h5::group g, std::string const &name, bravais_lattice &bl);
 
    private:
    matrix_t units_     = matrix_t::zeros(3, 3);
    matrix_t units_inv_ = matrix_t::zeros(3, 3);
    std::vector<r_t> atom_orb_pos_;
    std::vector<std::string> atom_orb_name_;
    int ndim_;
  };

 
} // namespace triqs::lattice