molecule.h

#pragma once

#include "auxiliary.h"
#include "particle.h"
#include "random.h"
#include <set>

namespace Faunus {

namespace pairpotential {
struct BondData;
}

namespace Geometry {
struct GeometryBase;
}

class QuaternionRotate;
class MoleculeData;

struct MoleculeInserter
{
    virtual ParticleVector operator()(const Geometry::GeometryBase& geo, MoleculeData& mol,
                                      const ParticleVector& other_particles) = 0;
    virtual void from_json(const json& j);
    virtual void to_json(json& j) const;
    virtual ~MoleculeInserter() = default;
};

void from_json(const json& j, MoleculeInserter& inserter);
void to_json(json& j, const MoleculeInserter& inserter);

class RandomInserter : public MoleculeInserter
{
  private:
    void translateRotateAtomicGroup(const Geometry::GeometryBase& geo,
                                    Faunus::QuaternionRotate& rotator,
                                    ParticleVector& particles) const;
    void translateRotateMolecularGroup(const Geometry::GeometryBase& geo, QuaternionRotate& rotator,
                                       ParticleVector& particles) const;

  public:
    Point dir = {1, 1, 1};    
    Point offset = {0, 0, 0}; 
    bool rotate = true; 
    bool keep_positions = false; 
    bool allow_overlap = false;  
    int max_trials = 20'000;     //!< Maximum number of container overlap checks

    ParticleVector
    operator()(const Geometry::GeometryBase& geo, MoleculeData& molecule,
               const ParticleVector& ignored_other_particles = ParticleVector()) override;
    void from_json(const json& j) override;
    void to_json(json& j) const override;
};

struct Conformation
{
    std::vector<Point> positions;
    std::vector<double> charges;
    bool empty() const;
    void copyTo(ParticleVector& particles) const; //!< Copy conformation into particle vector
};

class ExclusionsSimple
{
    bool any_exclusions = false; //!< true if at least one excluded interaction is present
    int size;                    //!< number of particles within the group; matrix = size × size
    using char_bool =
        unsigned char; //!< obs: std::vector<bool> uses packing with performance implications
    std::shared_ptr<std::vector<char_bool>>
        excluded_pairs; //!< 1D exclusion matrix; shared ptr saves copying

  public:
    using AtomPair = std::pair<int, int>;
    static ExclusionsSimple create(int atoms_cnt, const std::vector<std::pair<int, int>>& pairs);
    explicit ExclusionsSimple(int size = 0);
    void add(int i, int j);
    void add(const std::vector<AtomPair>& pairs);
    bool isExcluded(int i, int j) const;
    bool empty() const; //!< true if no excluded interactions at all
    friend void from_json(const json& j, ExclusionsSimple& exclusions);
    friend void to_json(json& j, const ExclusionsSimple& exclusions);
};

inline bool ExclusionsSimple::isExcluded(int i, int j) const
{
    if (i > j) {
        std::swap(i, j);
    }
    // use the bracket syntax to skip checks to speed-up reading
    return any_exclusions && static_cast<bool>((*excluded_pairs)[i * size + j]);
}

inline bool ExclusionsSimple::empty() const
{
    return !any_exclusions;
}

void from_json(const json& j, ExclusionsSimple& exclusions);
void to_json(json& j, const ExclusionsSimple& exclusions);

class ExclusionsVicinity
{
  public:
    using char_bool =
        unsigned char; //!< obs: std::vector<bool> uses packing with performance implications
    using AtomPair = std::pair<int, int>;

  private:
    int atoms_cnt = 0;         //!< count of atoms in the molecule; unmutable
    int max_bond_distance = 0; //!< max distance (difference) between indices of excluded particles
    std::shared_ptr<std::vector<char_bool>>
        excluded_pairs; //!< 1D exclusion matrix; shared ptr saves copying

    int toIndex(int i, int j) const;
    AtomPair fromIndex(int n) const;

  public:
    static ExclusionsVicinity create(int atoms_cnt, const std::vector<AtomPair>& pairs);
    explicit ExclusionsVicinity(int atoms_cnt = 0, int max_difference = 0);

    void add(int i, int j);
    void add(const std::vector<std::pair<int, int>>& pairs);
    bool isExcluded(int i,
                    int j)
        const; //!< @param i, j indices of atoms within molecule with excluded nonbonded interaction
    bool empty() const; //!< true if no excluded interactions at all
    friend void to_json(json& j, const ExclusionsVicinity& exclusions);
};

inline bool ExclusionsVicinity::isExcluded(int i, int j) const
{
    if (i > j) {
        std::swap(i, j);
    }
    // use bracket syntax to skip checks to speed-up reading
    return (j - i <= max_bond_distance && static_cast<bool>((*excluded_pairs)[toIndex(i, j)]));
}

inline bool ExclusionsVicinity::empty() const
{
    return max_bond_distance == 0;
}

inline int ExclusionsVicinity::toIndex(int i, int j) const
{
    return i * max_bond_distance + (j - i - 1);
}

inline ExclusionsVicinity::AtomPair ExclusionsVicinity::fromIndex(const int n) const
{
    auto i = n / max_bond_distance;
    auto j = n % max_bond_distance + i + 1;
    return {i, j};
}

// void from_json(const json &j, ExclusionsVicinity &exclusions);  // not implemented
void to_json(json& j, const ExclusionsVicinity& exclusions);

class MoleculeData
{
  public:
    using index_type = int;

  private:
    json json_cfg; //!< data useful only for to_json
    index_type _id = 0;
    bool implicit = false; //!< Is molecule implicit and explicitly absent from simulation cell?

  protected:
    ExclusionsVicinity exclusions; //!< Implementation of isPairExcluded;
  public:
    std::shared_ptr<MoleculeInserter> inserter = nullptr; //!< Functor for insertion into space

    index_type& id();                                 //!< Type id
    const index_type& id() const;                     //!< Type id
    void createMolecularConformations(const json& j); //!< Add conformations if appropriate
    void setConformationWeights(const json& j);       //!< Add weights for conformations

    std::string name;    //!< Molecule name
    bool atomic = false; //!< True if atomic group (salt etc.)
    bool compressible =
        false;             //!< True if compressible group (scales internally upon volume change)
    bool rigid = false;    //!< True if particle should be considered as rigid
    double activity = 0.0; //!< Chemical activity (mol/l)

    std::vector<AtomData::index_type> atoms; //!< Sequence of atoms in molecule (atom id's)
    BasePointerVector<pairpotential::BondData> bonds;
    WeightedDistribution<ParticleVector> conformations; 
    size_t numConformations() const;                    

    MoleculeData();
    MoleculeData(const std::string& name, const ParticleVector& particles,
                 const BasePointerVector<pairpotential::BondData>& bonds);

    bool isImplicit() const; 
    bool isPairExcluded(int i, int j) const;
    bool isMolecular() const;
    bool isAtomic() const;

    void setInserter(std::shared_ptr<MoleculeInserter> ins);

    ParticleVector getRandomConformation(Geometry::GeometryBase& geo,
                                         const ParticleVector& otherparticles = ParticleVector());

    friend class MoleculeBuilder;
    friend void to_json(json& j, const MoleculeData& a);
    friend void from_json(const json& j, MoleculeData& a);
}; // end of class

inline bool MoleculeData::isPairExcluded(int i, int j) const
{
    return exclusions.isExcluded(i, j);
}

void to_json(json& j, const MoleculeData& a);
void from_json(const json& j, MoleculeData& a);
void from_json(const json& j, std::vector<MoleculeData>& v);

// global instance of molecule vector
extern std::vector<MoleculeData> molecules;

struct UnknownMoleculeError : public GenericError
{
    explicit UnknownMoleculeError(std::string_view molecule_name);
};

MoleculeData& findMoleculeByName(std::string_view name);

std::vector<MoleculeData::index_type> parseMolecules(const json& j);

class MoleculeBuilder
{
    bool is_used = false;      
    std::string molecule_name; 
    ParticleVector particles;  
    decltype(MoleculeData::bonds) bonds;
    std::vector<std::pair<int, int>> exclusion_pairs;

  protected:
    static bool
    isFasta(const json& j); 
    void readCompoundValues(
        const json& j); 
    void readAtomic(const json& j); 
    void
    readParticles(const json& j); 
    void readBonds(const json& j);      
    void readFastaBonds(const json& j); 
    void readExclusions(const json& j); 
    static std::shared_ptr<MoleculeInserter> createInserter(const json& j);

  public:
    void from_json(const json& j, MoleculeData& molecule);
};

class MoleculeStructureReader
{
    bool read_charges; 
  protected:
    static void readArray(ParticleVector& particles, const json& j_particles);
    static void readFasta(ParticleVector& particles, const json& j);

  public:
    MoleculeStructureReader(bool read_charges = true);
    void readFile(ParticleVector& particles, const std::string& filename) const;
    void readJson(ParticleVector& particles, const json& j) const;
};

class NeighboursGenerator
{
  public:
    using AtomList = std::vector<int>;
    using AtomPair = std::pair<int, int>;
    using AtomPairList = std::vector<AtomPair>;
    using BondVector = decltype(MoleculeData::bonds);

  private:
    std::map<int, AtomList>
        bond_map; 
    std::vector<std::set<AtomList>> paths;

    void createBondMap(const BondVector& bonds); 
    void generatePaths(int bonded_distance);     

  public:
    NeighboursGenerator(const BondVector& bonds);
    void generatePairs(AtomPairList& pairs, int bond_distance);
};

class ReactionData
{
  public:
    using StoichiometryMap = std::map<int, int>; 
    using AtomicAndMolecularPair = std::pair<const StoichiometryMap&, const StoichiometryMap&>;
    using StoichiometryPair =
        StoichiometryMap::value_type; 
    using MapFilter = std::function<bool(const ReactionData::StoichiometryPair&)>;
    enum class Direction : char
    {
        LEFT = 0,
        RIGHT = 1
    };

    const static MapFilter is_implicit_group;
    const static MapFilter not_implicit_group;
    const static MapFilter not_implicit_atom;
    const static MapFilter is_atomic_group;
    const static MapFilter is_molecular_group;

  private:
    friend void from_json(const json&, ReactionData&);
    friend void to_json(json&, const ReactionData&);

    Direction direction = Direction::RIGHT;           
    std::vector<std::string> left_names, right_names; 
    StoichiometryMap left_molecules;                  
    StoichiometryMap right_molecules;                 
    StoichiometryMap left_atoms;                      
    StoichiometryMap right_atoms;                     
    double lnK = 0.0;            
    double lnK_unmodified = 0.0; 
    std::string reaction_str;    

    static std::pair<decltype(Faunus::atoms)::const_iterator,
                     decltype(Faunus::molecules)::const_iterator>
    findAtomOrMolecule(const std::string& atom_or_molecule_name);

  public:
    void setDirection(Direction);                
    void setRandomDirection(Random& random);     
    Direction getDirection() const;              
    void reverseDirection();                     
    AtomicAndMolecularPair getProducts() const;  
    AtomicAndMolecularPair getReactants() const; 
    double freeEnergy() const;       
    bool containsAtomicSwap() const; 
    const std::string& getReactionString() const; 

    std::pair<std::set<int>, std::set<int>>
    participatingAtomsAndMolecules() const; 

    bool only_neutral_molecules = false; 
}; 

std::pair<std::vector<std::string>, std::vector<std::string>>
parseReactionString(const std::string& process_string);

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

extern std::vector<ReactionData> reactions; // global instance

} // namespace Faunus