structures.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 <map>
#include <array>
#include <vector>
#include <utility>
#include <numeric>
#include <limits>
#include <boost/math/special_functions/pow.hpp>
#include <boost/mpi.hpp>
#include <Eigen/Dense>
#include <Eigen/StdVector>
#include <cmakevars.hpp>
#include <constants.hpp>
#include <utilities.hpp>
#include <exceptions.hpp>
#include <periodic_table.hpp>

namespace alma {
// Forward declarations placed here for convenience.
class Symmetry_operations;
class Gamma_grid;
class Threeph_process;

using Triple_int = std::array<int, 3>;

template <typename T>
using Triple_int_map =
    std::map<Triple_int, T, Container_comparator<Triple_int>>;

class Crystal_structure {
public:
    const Eigen::Matrix3d lattvec;
    const Eigen::Matrix<double, 3, Eigen::Dynamic> positions;
    const std::vector<std::string> elements;
    const double V;
    const Eigen::MatrixXd rlattvec;
    const std::vector<int> numbers;
    Crystal_structure(Eigen::Matrix3d _lattvec,
                      Eigen::Matrix<double, 3, Eigen::Dynamic> _positions,
                      std::vector<std::string> _elements,
                      std::vector<int> _numbers)
        : lattvec(std::move(_lattvec)), positions(std::move(_positions)),
          elements(std::move(_elements)),
          V(std::fabs(this->lattvec.determinant())),
          rlattvec(2. * constants::pi * this->lattvec.inverse().transpose()),
          numbers(std::move(_numbers)) {
        std::partial_sum(this->numbers.begin(),
                         this->numbers.end(),
                         std::back_inserter(this->partsums));

        for (int i = 0; i < this->get_natoms(); ++i) {
            this->masses.emplace_back(alma::get_mass(this->get_element(i)));
            this->gfactors.emplace_back(
                alma::get_gfactor(this->get_element(i)));
        }
    }

    Crystal_structure(Eigen::Matrix3d _lattvec,
                      Eigen::Matrix<double, 3, Eigen::Dynamic> _positions,
                      std::vector<std::string> _elements,
                      std::vector<int> _numbers,
                      std::vector<double> _masses)
        : Crystal_structure(_lattvec, _positions, _elements, _numbers) {
        if (_masses.size() != static_cast<std::size_t>(this->get_natoms())) {
            throw value_error("one mass per atom must be provided");
        }
        for (int i = 0; i < this->get_natoms(); ++i) {
            if (_masses[i] <= 0.) {
                throw value_error("all masses must be positive");
            }
            this->masses[i] = _masses[i];
        }
    }

    inline std::vector<std::string>::size_type get_nelements() const {
        return this->elements.size();
    }


    inline int get_natoms() const {
        return this->positions.cols();
    }
    inline std::string get_element(int i) const {
        if (i >= this->get_natoms())
            throw value_error("invalid atom number");
        return this->elements[std::distance(
            this->partsums.begin(),
            std::lower_bound(
                this->partsums.begin(), this->partsums.end(), i + 1))];
    }


    inline double get_mass(int i) const {
        if (i >= this->get_natoms())
            throw value_error("invalid atom number");
        return this->masses[i];
    }


    inline double get_gfactor(int i) const {
        if (i >= this->get_natoms())
            throw value_error("invalid atom number");
        return this->gfactors[i];
    }


    inline bool is_alloy() const {
        for (auto& e : this->elements)
            if (e.find(";") != std::string::npos)
                return true;

        return false;
    }

    inline Eigen::MatrixXd map_to_firstbz(
        const Eigen::Ref<const Eigen::Vector3d>& q) const {
        constexpr int sbound = 1;
        // Step 1: find a single image using a standard approach
        Eigen::Vector3d qbz{this->rlattvec.colPivHouseholderQr().solve(q)};
        qbz -= qbz.array().round().matrix().eval();
        qbz = (this->rlattvec * qbz).eval();
        double n2 = qbz.squaredNorm();
        for (int i = -sbound; i <= sbound; ++i) {
            for (int j = -sbound; j <= sbound; ++j) {
                for (int k = -sbound; k <= sbound; ++k) {
                    Eigen::Vector3d qbzp{qbz + i * this->rlattvec.col(0) +
                                         j * this->rlattvec.col(1) +
                                         k * this->rlattvec.col(2)};
                    double n2p = qbzp.squaredNorm();
                    if (n2p < n2) {
                        n2 = n2p;
                        qbz = qbzp;
                    }
                }
            }
        }
        // Step 2: perform a search to find equivalent images.
        std::size_t found = 0;
        int ncols = boost::math::pow<3>(2 * sbound + 1);
        Eigen::MatrixXd nruter(3, ncols);
        for (int i = 0; i <= sbound; ++i) {
            for (int j = -sbound; j <= sbound; ++j) {
                for (int k = -sbound; k <= sbound; ++k) {
                    Eigen::Vector3d qbzp{qbz + i * this->rlattvec.col(0) +
                                         j * this->rlattvec.col(1) +
                                         k * this->rlattvec.col(2)};
                    double n2p = qbzp.squaredNorm();
                    if (alma::almost_equal(n2, n2p)) {
                        nruter.col(found) = qbzp;
                        ++found;
                    }
                }
            }
        }
        return nruter.leftCols(found);
    }


private:
    std::vector<int> partsums;
    std::vector<double> masses;
    std::vector<double> gfactors;
};

template <class T> class General_harmonic_ifcs {
public:
    const std::vector<Triple_int> pos;
    const std::vector<T> ifcs;
    const int na;
    const int nb;
    const int nc;
    General_harmonic_ifcs(std::vector<Triple_int> _pos,
                          std::vector<T> _ifcs,
                          int _na,
                          int _nb,
                          int _nc)
        : pos(std::move(_pos)), ifcs(std::move(_ifcs)), na(_na), nb(_nb),
          nc(_nc) {
    }


    // Return the number of unit cells.
    inline std::vector<std::string>::size_type get_ncells() const {
        return this->pos.size();
    }
};

using Harmonic_ifcs = General_harmonic_ifcs<Eigen::MatrixXd>;

class Dielectric_parameters {
public:
    const std::vector<Eigen::MatrixXd> born;
    const Eigen::Matrix3d epsilon;
    Dielectric_parameters(std::vector<Eigen::MatrixXd> _born,
                          Eigen::Matrix3d _epsilon)
        : born(std::move(_born)), epsilon(std::move(_epsilon)) {
    }


    Dielectric_parameters() {
    }
};

class Supercell_index {
public:
    const int index;
    const int ia;
    const int ib;
    const int ic;
    const int iatom;
    // Basic constructor.
    Supercell_index(const int _index,
                    const int _ia,
                    const int _ib,
                    const int _ic,
                    const int _iatom)
        : index(_index), ia(_ia), ib(_ib), ic(_ic), iatom(_iatom) {
    }


    inline const Triple_int get_pos() const {
        return Triple_int({{this->ia, this->ib, this->ic}});
    }
};

class Supercell_index_builder {
public:
    const int na;
    const int nb;
    const int nc;
    const int natoms;
    Supercell_index_builder(const int _na,
                            const int _nb,
                            const int _nc,
                            const int _natoms);
    Supercell_index create_index(const int index) const;

    Supercell_index create_index(const int ia,
                                 const int ib,
                                 const int ic,
                                 const int iatom) const;

    Supercell_index create_index_safely(const int index) const;

    Supercell_index create_index_safely(const int ia,
                                        const int ib,
                                        const int ic,
                                        const int iatom) const;
};

class Thirdorder_ifcs {
public:
    const Eigen::VectorXd rj;
    const Eigen::VectorXd rk;
    const std::size_t i;
    const std::size_t j;
    const std::size_t k;
    double& ifc(std::size_t alpha, std::size_t beta, std::size_t gamma) {
        return this->ifcs[gamma + 3 * (beta + 3 * alpha)];
    }


    const double& ifc(std::size_t alpha,
                      std::size_t beta,
                      std::size_t gamma) const {
        return this->ifcs[gamma + 3 * (beta + 3 * alpha)];
    }


    Thirdorder_ifcs(const Eigen::VectorXd& _rj,
                    const Eigen::VectorXd& _rk,
                    std::size_t _i,
                    std::size_t _j,
                    std::size_t _k)
        : rj(std::move(_rj)), rk(std::move(_rk)), i(_i), j(_j), k(_k) {
    }


private:
    std::array<double, 27> ifcs;
};
} // namespace alma