qpoint_grid.hpp

Go to the documentation of this file.

// Copyright 2015-2018 The ALMA Project Developers
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//   http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
// implied. See the License for the specific language governing
// permissions and limitations under the License.

#pragma once


#include <boost/mpi.hpp>
#include <structures.hpp>
#include <dynamical_matrix.hpp>

namespace alma {
using Tetrahedron = std::array<std::size_t, 4>;
using Triangle = std::array<std::size_t, 3>;

class Gamma_grid {
public:
    const int na;
    const int nb;
    const int nc;
    const std::size_t nqpoints;
    const Eigen::MatrixXd rlattvec;
    const Eigen::MatrixXd dq;
    Gamma_grid(const Crystal_structure& poscar,
               const Symmetry_operations& symms,
               const Harmonic_ifcs& force_constants,
               int _na,
               int _nb,
               int _nc);


    Gamma_grid(const Crystal_structure& poscar,
               const Symmetry_operations& symms,
               const Harmonic_ifcs& force_constants,
               const Dielectric_parameters& born,
               int _na,
               int _nb,
               int _nc);


    void enforce_asr() {
        this->spectrum[0].omega.head(3).fill(0.);
    }


    std::size_t three_to_one(const std::array<int, 3>& indices) const {
        auto ia = python_mod(indices[0], this->na);
        auto ib = python_mod(indices[1], this->nb);
        auto ic = python_mod(indices[2], this->nc);

        return ic + nc * (ib + nb * ia);
    }


    std::array<int, 3> one_to_three(std::size_t iq) const {
        std::array<int, 3> nruter;
        nruter[2] = iq % nc;
        iq /= nc;
        nruter[1] = iq % nb;
        nruter[0] = iq / nb;
        return nruter;
    }


    std::size_t get_nequivalences() const {
        return this->equivalences.size();
    }


    const Spectrum_at_point& get_spectrum_at_q(int iq) const {
        std::size_t index = python_mod(iq, this->nqpoints);

        return this->spectrum[index];
    }


    std::size_t get_cardinal(std::size_t ic) const {
        if (ic > this->equivalences.size())
            throw value_error("wrong equivalence class index");
        return this->equivalences[ic].size();
    }


    std::size_t get_representative(std::size_t ic) const {
        if (ic > this->equivalences.size())
            throw value_error("wrong equivalence class index");
        return this->equivalences[ic][0];
    }


    std::vector<std::size_t> get_equivalence(std::size_t ic) const {
        if (ic > this->equivalences.size())
            throw value_error("wrong equivalence class index");
        return this->equivalences[ic];
    }


    Eigen::VectorXd get_q(std::size_t iq) const {
        std::size_t index = python_mod(iq, this->nqpoints);

        return this->cpos.col(index);
    }


    double base_sigma(const Eigen::VectorXd& v) const {
        return (v.transpose() * this->dq).norm() / std::sqrt(12.);
    }


    Gamma_grid(const Crystal_structure& poscar,
               const Symmetry_operations& symms,
               int _na,
               int _nb,
               int _nc);


    std::size_t polar_opposite(std::size_t q) {
        auto indices = this->one_to_three(q);

        return this->three_to_one({{-indices[0], -indices[1], -indices[2]}});
    }


    std::vector<size_t> equivalent_qpoints(std::size_t original) const {
        std::size_t index = python_mod(original, this->nqpoints);
        return this->symmetry_map[index];
    }


    std::vector<std::array<std::size_t, 2>> equivalent_qpairs(
        const std::array<std::size_t, 2>& original) const;


    std::vector<std::array<std::size_t, 3>> equivalent_qtriplets(
        const std::array<std::size_t, 3>& original) const;


    std::vector<Tetrahedron> get_tetrahedra(std::size_t q) const {
        // Find out the indices of the eight corners of the
        // microcell.
        std::size_t i000 = q;
        auto indices = this->one_to_three(i000);
        std::size_t i001 =
            this->three_to_one({{indices[0], indices[1], indices[2] + 1}});
        std::size_t i010 =
            this->three_to_one({{indices[0], indices[1] + 1, indices[2]}});
        std::size_t i011 =
            this->three_to_one({{indices[0], indices[1] + 1, indices[2] + 1}});
        std::size_t i100 =
            this->three_to_one({{indices[0] + 1, indices[1], indices[2]}});
        std::size_t i101 =
            this->three_to_one({{indices[0] + 1, indices[1], indices[2] + 1}});
        std::size_t i110 =
            this->three_to_one({{indices[0] + 1, indices[1] + 1, indices[2]}});
        std::size_t i111 = this->three_to_one(
            {{indices[0] + 1, indices[1] + 1, indices[2] + 1}});

        // Build the four equivalent tetrahedra.
        std::vector<Tetrahedron> nruter;
        nruter.reserve(5);
        nruter.emplace_back(Tetrahedron({{i000, i001, i010, i100}}));
        nruter.emplace_back(Tetrahedron({{i110, i111, i100, i010}}));
        nruter.emplace_back(Tetrahedron({{i101, i100, i111, i001}}));
        nruter.emplace_back(Tetrahedron({{i011, i010, i001, i111}}));
        // And the central, inequivalent one.
        nruter.emplace_back(Tetrahedron({{i010, i100, i111, i001}}));
        return nruter;
    }


    std::vector<Triangle> get_triangles(std::size_t ia) const;


    std::size_t getParentIdx(std::size_t iq) const;


    std::size_t getSymIdxToParent(std::size_t iq) const;


    std::vector<std::vector<std::size_t>> getSymmetryMap() {
        return this->symmetry_map;
    }


    template <typename T>
    auto copy_symmetry(std::size_t iq,
                       const Symmetry_operations& symms,
                       const Eigen::MatrixBase<T>& x) const
        -> Eigen::Matrix<typename T::Scalar, Eigen::Dynamic, Eigen::Dynamic> {
        std::size_t index = python_mod(iq, this->nqpoints);
        std::size_t nops = symms.get_nsym();

        typename T::PlainObject nruter(x);
        nruter.fill(0.);
        std::size_t nfound = 0;

        for (std::size_t iop = 0; iop < nops; ++iop) {
            if (this->symmetry_map[index][2 * iop] == index) {
                nruter += symms.rotate_v(x, iop, true);
                ++nfound;
            }
        }
        nruter /= nfound;
        return nruter;
    }

private:
    friend void save_bulk_hdf5(const char* filename,
                               const std::string& description,
                               const Crystal_structure& cell,
                               const Symmetry_operations& symmetries,
                               const Gamma_grid& grid,
                               const std::vector<Threeph_process>& processes,
                               const boost::mpi::communicator& comm);

    friend std::tuple<std::string,
                      std::unique_ptr<Crystal_structure>,
                      std::unique_ptr<Symmetry_operations>,
                      std::unique_ptr<Gamma_grid>,
                      std::unique_ptr<std::vector<Threeph_process>>>
    load_bulk_hdf5(const char* filename, const boost::mpi::communicator& comm);

    Eigen::MatrixXd cpos;
    std::vector<Spectrum_at_point> spectrum;
    std::vector<std::vector<std::size_t>> equivalences;
    std::vector<std::vector<std::size_t>> symmetry_map;
    void initialize_cpos();

    std::vector<Spectrum_at_point> compute_my_spectrum(
        const Dynamical_matrix_builder& factory,
        const Symmetry_operations& symms,
        const boost::mpi::communicator& communicator);

    void fill_equivalences(const Symmetry_operations& symms);

    std::vector<std::size_t> parentlookup;

    void fill_parentlookup();

    void fill_map(const Symmetry_operations& symms);
};
} // namespace alma