speciation.h

#pragma once

#include "move.h"

namespace Faunus::Speciation {

class ReactionValidator
{
  private:
    const Space& spc;
    [[nodiscard]] bool canSwapAtoms(const ReactionData& reaction) const;
    [[nodiscard]] bool canReduceImplicitGroups(const ReactionData& reaction) const;
    [[nodiscard]] bool canReduceMolecularGroups(const ReactionData& reaction) const;
    [[nodiscard]] bool canProduceMolecularGroups(const ReactionData& reaction) const;
    [[nodiscard]] bool canReduceAtomicGroups(const ReactionData& reaction) const;
    [[nodiscard]] bool canProduceAtomicGroups(const ReactionData& reaction) const;

  public:
    explicit ReactionValidator(const Space& spc);
    bool isPossible(const ReactionData& reaction) const; 
};

class GroupDeActivator
{
  public:
    using ChangeAndBias = std::pair<Change::GroupChange, double>; 
    using OptionalInt = std::optional<int>;
    virtual ChangeAndBias activate(Group& group, OptionalInt num_particles = std::nullopt) = 0;
    virtual ChangeAndBias deactivate(Group& group, OptionalInt num_particles = std::nullopt) = 0;
    virtual ~GroupDeActivator() = default;
};

class AtomicGroupDeActivator : public GroupDeActivator
{
  private:
    Space& spc;     
    Space& old_spc; 
    Random& random;

  public:
    AtomicGroupDeActivator(Space& spc, Space& old_spc, Random& random);
    ChangeAndBias activate(Group& group, OptionalInt number_to_insert) override;
    ChangeAndBias deactivate(Group& group, OptionalInt number_to_delete) override;
};

class MolecularGroupDeActivator : public GroupDeActivator
{
  private:
    Space& spc;
    Random& random;
    const bool apply_bond_bias; 
    [[nodiscard]] double getBondEnergy(const Group& group) const;
    virtual void
    setPositionAndOrientation(Group& group) const; 

  public:
    MolecularGroupDeActivator(Space& spc, Random& random, bool apply_bond_bias);
    ChangeAndBias activate(Group& group, OptionalInt num_particles = std::nullopt) override;
    ChangeAndBias deactivate(Group& group, OptionalInt num_particles = std::nullopt) override;
};

class ReactionDirectionRatio
{
  private:
    using reaction_iterator = decltype(Faunus::reactions)::iterator;

    struct AcceptanceData
    {
        Average<double> right; 
        Average<double> left;  
        void update(ReactionData::Direction direction, bool accept);
    };

    std::map<reaction_iterator, AcceptanceData> acceptance;

  public:
    AcceptanceData& operator[](reaction_iterator iter);
    void to_json(json& j) const;
};

} // end of namespace Faunus::Speciation

namespace Faunus::move {

class SpeciationMove : public Move
{
  private:
    using reaction_iterator = decltype(Faunus::reactions)::iterator;
    Random random_internal;     
    reaction_iterator reaction; 
    double bias_energy = 0.0;   
    Speciation::ReactionValidator reaction_validator;   
    Speciation::ReactionDirectionRatio direction_ratio; 
    std::unique_ptr<Speciation::GroupDeActivator>
        molecular_group_bouncer; 
    std::unique_ptr<Speciation::GroupDeActivator>
        atomic_group_bouncer; 

    std::map<MoleculeData::index_type, Average<double>>
        average_reservoir_size; 

    void _to_json(json& j) const override;
    void _from_json(const json& j) override;
    void _move(Change& change) override;
    void _accept(Change& change) override;
    void _reject(Change& change) override;

    void setRandomReactionAndDirection();     
    void atomicSwap(Change& change);          
    void deactivateReactants(Change& change); 
    void activateProducts(Change& change);    
    void deactivateAtomicGroups(Change& change);
    void activateAtomicGroups(Change& change);
    void deactivateMolecularGroups(Change& change);
    void activateMolecularGroups(Change& change);
    void
    updateGroupMassCenters(const Change& change) const; 
    SpeciationMove(Space& spc, Space& old_spc, std::string_view name, std::string_view cite);

  public:
    SpeciationMove(Space& spc, Space& old_spc);
    double bias(Change& change, double old_energy, double new_energy) override;
    static void swapParticleProperties(Particle& particle, int new_atomid);
};

} // namespace Faunus::move