processes.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 <cstddef>
#include <cmath>
#include <array>
#include <vector>
#include <functional>
#include <boost/mpi.hpp>
#include <boost/math/distributions/normal.hpp>
#include <boost/serialization/array.hpp>
#include <structures.hpp>
#include <qpoint_grid.hpp>

// Forward declarations of elements documented later on.
// Note that serialization of std::vectors of objects without
// a default constructor is broken in Boost version 1.58.0. See:
// http://stackoverflow.com/a/30437359/85371
// https://svn.boost.org/trac/boost/ticket/11342
// This will prevent ALMA from compiling. We check for that
// specific version in our cmake configuration.
namespace alma {
class Threeph_process;
}
namespace boost {
namespace serialization {
template <class Archive>
inline void save_construct_data(Archive& ar,
                                const alma::Threeph_process* t,
                                const unsigned int file_version);
}
} // namespace boost

namespace alma {
enum class threeph_type { emission = -1, absorption = 1 };
class Threeph_process {
private:
    friend class boost::serialization::access;
    template <class Archive>
    friend void boost::serialization::save_construct_data(
        Archive& ar,
        const Threeph_process* t,
        const unsigned int file_version);

    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);

    const double domega;
    const double sigma;
    bool gaussian_computed;
    double gaussian;
    bool vp2_computed;
    double vp2;
    template <class Archive>
    void serialize(Archive& ar, const unsigned int file_version) {
        ar & this->gaussian_computed;

        if (this->gaussian_computed)
            ar & this->gaussian;
        ar & this->vp2_computed;

        if (this->vp2_computed)
            ar & this->vp2;
    }


public:
    const std::size_t c;
    const std::array<std::size_t, 3> q;
    const std::array<std::size_t, 3> alpha;
    const threeph_type type;
    Threeph_process(std::size_t _c,
                    const std::array<std::size_t, 3>& _q,
                    const std::array<std::size_t, 3>& _alpha,
                    threeph_type _type,
                    double _domega,
                    double _sigma)
        : domega(_domega), sigma(_sigma), gaussian_computed(false),
          vp2_computed(false), c(_c), q(std::move(_q)),
          alpha(std::move(_alpha)), type(_type) {
    }


    Threeph_process(const Threeph_process& original)
        : domega(original.domega), sigma(original.sigma),
          gaussian_computed(original.gaussian_computed),
          gaussian(original.gaussian), vp2_computed(original.vp2_computed),
          vp2(original.vp2), c(original.c), q(original.q),
          alpha(original.alpha), type(original.type) {
    }


    double compute_gaussian() {
        if (!this->gaussian_computed) {
            auto distr = boost::math::normal(0., this->sigma);
            this->gaussian = boost::math::pdf(distr, this->domega);
            this->gaussian_computed = true;
        }
        return this->gaussian;
    }


    double compute_weighted_gaussian(const Gamma_grid& grid, double T);


    double compute_vp2(const Crystal_structure& cell,
                       const Gamma_grid& grid,
                       const std::vector<Thirdorder_ifcs>& thirdorder);

    inline double get_vp2() const {
        if (!this->vp2_computed)
            throw exception("vp2 is not available");
        return this->vp2;
    }


    inline bool is_vp2_computed() const {
        return this->vp2_computed;
    }


    double compute_gamma(const Gamma_grid& grid, double T);

    double compute_gamma_reduced(const Gamma_grid& grid, double T);

    Eigen::ArrayXd compute_collision(const Gamma_grid& grid,
                                     const Eigen::ArrayXXd& n0);
};
std::vector<Threeph_process> find_allowed_threeph(
    const Gamma_grid& grid,
    const boost::mpi::communicator& communicator,
    double scalebroad = 1.0);

Eigen::ArrayXXd calc_w0_threeph(const alma::Gamma_grid& grid,
                                std::vector<alma::Threeph_process>& processes,
                                double T,
                                const boost::mpi::communicator& comm);

Eigen::ArrayXXd calc_w0_threeph(
    const alma::Gamma_grid& grid,
    std::vector<alma::Threeph_process>& processes,
    double T,
    std::function<bool(const Threeph_process&)> filter,
    const boost::mpi::communicator& comm);
} // namespace alma

namespace boost {
namespace serialization {
template <class Archive>
inline void save_construct_data(Archive& ar,
                                const alma::Threeph_process* t,
                                const unsigned int file_version) {
    ar << t->c;
    ar << make_array(t->q.data(), t->q.size());
    ar << make_array(t->alpha.data(), t->alpha.size());
    ar << t->type;
    ar << t->domega;
    ar << t->sigma;
}


template <class Archive>
inline void load_construct_data(Archive& ar,
                                alma::Threeph_process* t,
                                const unsigned int file_version) {
    std::size_t c;

    ar >> c;
    std::array<std::size_t, 3> q;
    ar >> make_array(q.data(), q.size());
    std::array<std::size_t, 3> alpha;
    ar >> make_array(alpha.data(), alpha.size());
    alma::threeph_type type;
    ar >> type;
    double domega;
    ar >> domega;
    double sigma;
    ar >> sigma;
    ::new (t) alma::Threeph_process(c, q, alpha, type, domega, sigma);
}
} // namespace serialization
} // namespace boost