particle.h

#pragma once
#include "core.h"
#include "atomdata.h"
#include "tensor.h"
#include <spdlog/spdlog.h>

namespace Eigen {
using Matrix3d = Matrix<double, 3, 3>;
using Quaterniond = Quaternion<double>;
} // namespace Eigen

namespace Faunus {

struct ParticlePropertyBase
{
    virtual void to_json(json& j) const = 0;                          
    virtual void from_json(const json& j) = 0;                        
    void rotate(const Eigen::Quaterniond& q, const Eigen::Matrix3d&); 
    virtual ~ParticlePropertyBase() = default;
};

template <typename... Ts> auto to_json(json&) -> typename std::enable_if<sizeof...(Ts) == 0>::type
{
} 

template <typename T, typename... Ts>
void to_json(json& j, const ParticlePropertyBase& a, const Ts&... rest)
{
    a.to_json(j);
    to_json<Ts...>(j, rest...);
} 

// JSON --> Particle
template <typename... Ts>
auto from_json(const json&) -> typename std::enable_if<sizeof...(Ts) == 0>::type
{
}

template <typename T, typename... Ts>
void from_json(const json& j, ParticlePropertyBase& a, Ts&... rest)
{
    a.from_json(j);
    from_json<Ts...>(j, rest...);
}

struct Radius : public ParticlePropertyBase
{
    double radius = 0.0; 
    void to_json(json& j) const override;
    void from_json(const json& j) override;

    template <class Archive> void serialize(Archive& archive) { archive(radius); }

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
}; 

struct Charge : public ParticlePropertyBase
{
    double charge = 0.0; 
    void to_json(json& j) const override;
    void from_json(const json& j) override;

    template <class Archive> void serialize(Archive& archive) { archive(charge); }

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
}; 

struct Dipole : public ParticlePropertyBase
{
    Point mu = {0.0, 0.0, 0.0};                                       
    double mulen = 0.0;                                               
    void rotate(const Eigen::Quaterniond& q, const Eigen::Matrix3d&); 
    void to_json(json& j) const override;
    void from_json(const json& j) override;
    bool isDipolar() const;

    template <class Archive> void serialize(Archive& archive) { archive(mu, mulen); }

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};

struct Polarizable : public ParticlePropertyBase
{
    Point mui = {1.0, 0.0, 0.0};   
    double muilen = 0.0;           
    Tensor alpha = Tensor::Zero(); 
    void rotate(const Eigen::Quaterniond& q,
                const Eigen::Matrix3d& m); 
    void to_json(json& j) const override;
    void from_json(const json& j) override;
    bool isPolarizable() const; 

    template <class Archive> void serialize(Archive& archive) { archive(mui, muilen, alpha); }

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};

struct Quadrupole : public ParticlePropertyBase
{
    Tensor Q = Tensor::Zero(); 
    void rotate(const Eigen::Quaterniond& q,
                const Eigen::Matrix3d& m); 
    void to_json(json& j) const override;
    void from_json(const json& j) override;
    bool isQuadrupolar() const; 

    template <class Archive> void serialize(Archive& archive) { archive(Q); }

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
}; // Quadrupole property

class Cigar : public ParticlePropertyBase
{
  private:
  public:
    using string_t = nlohmann::basic_json<>::string_t;
    Point scdir = {1.0, 0.0, 0.0};    
    Point patchdir = {0.0, 1.0, 0.0}; 
    std::array<Point, 2> patchsides;
    double half_length = 0.0; 
    double pcanglsw = 0.0;    
    double pcangl = 0.0;      
    void rotate(const Eigen::Quaterniond& quaternion,
                const Eigen::Matrix3d& rotation_matrix); 
    void to_json(json& j) const override;
    void from_json(const json& j) override;
    void setDirections(
        const SpheroCylinderData& psc_data, const Point& new_direction,
        const Point& new_patch_direction); // initialize; run at start and after patch changes

    template <class Archive> void serialize(Archive& archive)
    {
        archive(scdir, patchdir, patchsides.at(0), patchsides.at(1));
    }

    bool isCylindrical() const; 

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
}; 

template <typename... Properties> class ParticleTemplate : public Properties...
{
  private:
    // rotate internal coordinates
    template <typename... Ts>
    auto _rotate(const Eigen::Quaterniond&, const Eigen::Matrix3d&) ->
        typename std::enable_if<sizeof...(Ts) == 0>::type
    {
    }

    template <typename T, typename... Ts>
    void _rotate(const Eigen::Quaterniond& q, const Eigen::Matrix3d& m, T& a, Ts&... rest)
    {
        a.rotate(q, m);
        _rotate<Ts...>(q, m, rest...);
    }

  public:
    ParticleTemplate()
        : Properties()... {};

    explicit ParticleTemplate(const AtomData& a)
        : Properties()...
    {
        *this = json(a).front();
    }

    void rotate(const Eigen::Quaterniond& q, const Eigen::Matrix3d& m)
    {
        _rotate<Properties...>(q, m, dynamic_cast<Properties&>(*this)...);
    } 

    template <class Archive> void serialize(Archive& archive)
    {
        __serialize<Properties...>(archive, dynamic_cast<Properties&>(*this)...);
    }

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};

template <typename... Properties> void to_json(json& j, const ParticleTemplate<Properties...>& a)
{
    to_json<Properties...>(j, Properties(a)...);
}

template <typename... Properties> void from_json(const json& j, ParticleTemplate<Properties...>& a)
{
    from_json<Properties...>(j, dynamic_cast<Properties&>(a)...);
}

class Particle
{
  public:
    using ParticleExtension = ParticleTemplate<Dipole, Quadrupole, Cigar>;
    std::unique_ptr<ParticleExtension> ext; 
    int id = -1;                            
    double charge = 0.0;                    
    Point pos = {0.0, 0.0, 0.0};            

    Particle() = default;
    Particle(const AtomData& a, const Point& pos);
    Particle(const AtomData& a);          
    Particle(const Particle&);            
    Particle& operator=(const Particle&); 
    const AtomData& traits() const;       
    void rotate(const Eigen::Quaterniond& quaternion,
                const Eigen::Matrix3d& rotation_matrix); 
    bool hasExtension() const;            
    ParticleExtension& createExtension(); 

    inline ParticleExtension& getExt()
    {
        return ext ? *ext : createExtension();
    } 

    inline const ParticleExtension& getExt() const
    {
        assert(ext);
        return *ext;
    } 
};

using ParticleVector = std::vector<Particle>;

template <class T>
concept RequireParticles =
    std::ranges::range<T> && std::is_convertible_v<std::ranges::range_value_t<T>, Particle>;

template <class T>
concept RequireParticleIterator = std::is_convertible_v<std::iter_value_t<T>, Particle>;

void from_json(const json&, Particle&);
void to_json(json&, const Particle&);

template <RequireParticleIterator Iterator>
void applyAtomDataCharges(Iterator first, Iterator last, const double max_difference = 1e-9)
{
    size_t mismatch_counter = 0;
    std::for_each(first, last, [&](Particle& particle) {
        const auto topology_charge = Faunus::atoms.at(particle.id).charge;
        if (std::fabs(topology_charge - particle.charge) > max_difference) {
            particle.charge = topology_charge;
            mismatch_counter++;
        }
    });
    if (mismatch_counter > 0) {
        faunus_logger->debug("{} charge(s) reset with topology values (from atomlist)",
                             mismatch_counter);
    }
}

} // namespace Faunus