symmetry.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 <structures.hpp>

namespace alma {
class Transformed_pair {
public:
    const std::array<std::size_t, 2> pair;
    const std::size_t operation;
    const bool exchange;
    Transformed_pair(std::size_t a,
                     std::size_t b,
                     std::size_t _operation,
                     bool _exchange)
        : pair({{a, b}}), operation(_operation), exchange(_exchange) {
    }
};


class Symmetry_operations {
public:
    const double symprec;
    Symmetry_operations(const Crystal_structure& structure,
                        double _symprec = 1e-5);
    std::size_t get_nsym() const {
        return this->translations.size();
    }
    int get_spacegroup_number() const {
        return this->sg_number;
    }
    std::string get_spacegroup_symbol() const {
        return this->sg_symbol;
    }


    std::string get_wyckoff() const {
        return this->wyckoff;
    }
    std::vector<int> get_equivalences() const {
        return this->equivalences;
    }


    template <typename T>
    auto transform_v(const Eigen::MatrixBase<T>& vector,
                     std::size_t index,
                     bool cartesian = false) const
        -> Eigen::Matrix<typename T::Scalar, Eigen::Dynamic, Eigen::Dynamic> {
        if (index >= this->translations.size())
            throw alma::value_error("invalid operation index");

        if (cartesian)
            return this->crotations[index] * vector +
                   this->ctranslations[index];
        else
            return this->rotations[index] * vector + this->translations[index];
    }


    template <typename T>
    auto rotate_v(const Eigen::MatrixBase<T>& vector,
                  std::size_t index,
                  bool cartesian = false) const
        -> Eigen::Matrix<typename T::Scalar, Eigen::Dynamic, Eigen::Dynamic> {
        if (index >= this->translations.size())
            throw alma::value_error("invalid operation index");

        if (cartesian)
            return this->crotations[index] * vector;
        else
            return this->rotations[index] * vector;
    }


    template <typename T>
    Eigen::Matrix<T, 3, 3> rotate_m(
        const Eigen::Ref<const Eigen::Matrix<T, 3, 3>>& matrix,
        std::size_t index,
        bool cartesian = false) const {
        if (index >= this->translations.size())
            throw alma::value_error("invalid operation index");

        if (cartesian)
            return this->crotations[index] * matrix *
                   this->crotations[index].transpose();
        else
            return this->rotations[index] * matrix *
                   this->rotations[index].transpose();
    }


    template <typename T>
    Eigen::Matrix<T, 3, 3> unrotate_m(
        const Eigen::Ref<const Eigen::Matrix<T, 3, 3>>& matrix,
        std::size_t index,
        bool cartesian = false) const {
        if (index >= this->translations.size())
            throw alma::value_error("invalid operation index");

        if (cartesian)
            return this->crotations[index].transpose() * matrix *
                   this->crotations[index];
        else
            return this->rotations[index].transpose() * matrix *
                   this->rotations[index];
    }


    std::size_t map_atom(std::size_t original, std::size_t index) const {
        if (original >= this->symmetry_map.size())
            throw alma::value_error("invalid atom index");

        if (index >= this->translations.size())
            throw alma::value_error("invalid operation index");
        return this->symmetry_map[original][index];
    }


    bool is_inversion(std::size_t index) const {
        return this->determinants[index] < 0.;
    }


    template <typename T>
    auto rotate_q(const Eigen::MatrixBase<T>& vector, std::size_t index) const
        -> Eigen::Matrix<typename T::Scalar, Eigen::Dynamic, Eigen::Dynamic> {
        if (index >= this->translations.size())
            throw alma::value_error("invalid operation index");
        return this->qrotations[index] * vector;
    }

    template <typename T>
    Eigen::Matrix<T, 3, 3> symmetrize_m(
        const Eigen::Ref<const Eigen::Matrix<T, 3, 3>>& matrix,
        bool cartesian = false) const {
        Eigen::Matrix<T, 3, 3> nruter;
        nruter.fill(0.);
        const std::vector<Eigen::MatrixXd>* rots;
        if (cartesian) {
            rots = &this->crotations;
        }
        else {
            rots = &this->rotations;
        }
        for (const auto& r : *rots) {
            nruter += r.transpose() * matrix * r;
        }
        nruter /= rots->size();
        return nruter;
    }

    std::vector<std::vector<Transformed_pair>> get_pair_classes() const;

private:
    std::vector<Eigen::VectorXd> translations;
    std::vector<Eigen::VectorXd> ctranslations;
    std::vector<Eigen::MatrixXd> rotations;
    std::vector<Eigen::MatrixXd> crotations;
    std::vector<Eigen::MatrixXd> qrotations;
    std::vector<double> determinants;
    std::vector<Eigen::MatrixXd> displacements;
    int sg_number;
    std::string sg_symbol;
    std::string wyckoff;
    std::vector<int> equivalences;
    std::vector<std::vector<std::size_t>> symmetry_map;
    void fill_map(const Crystal_structure& structure);
};
} // namespace alma