forcemove.h

#pragma once
#include "move.h"
#include "energy.h"

namespace Faunus::move {

class NormalRandomVector
{
    std::normal_distribution<double> normal_distribution;

  public:
    explicit NormalRandomVector(double mean = 0.0, double stddev = 1.0)
        : normal_distribution(mean, stddev) {};

    template <typename random_engine> Point operator()(random_engine& engine)
    {
        return {normal_distribution(engine), normal_distribution(engine),
                normal_distribution(engine)};
    }
};

class IntegratorBase
{
  protected:
    Space& spc;
    Energy::EnergyTerm& energy;
    IntegratorBase(Space&, Energy::EnergyTerm&);
    virtual ~IntegratorBase() = default;

  public:
    virtual void step(PointVector& velocities,
                      PointVector& forces) = 0; // todo shall we return real time progress?
    virtual void from_json(const json& j) = 0;
    virtual void to_json(json& j) const = 0;
};

void from_json(const json& j, IntegratorBase& i);
void to_json(json& j, const IntegratorBase& i);

class LangevinVelocityVerlet : public IntegratorBase
{
    double time_step;            
    double friction_coefficient; 
    NormalRandomVector
        random_vector; 
    [[nodiscard]] inline Point positionIncrement(
        const Point& velocity) const; 
    [[nodiscard]] inline Point
    velocityIncrement(const Point& force,
                      double mass) const; 
    inline Point velocityFluctuationDissipation(const Point& velocity, double mass);

  public:
    LangevinVelocityVerlet(Space& spc, Energy::EnergyTerm& energy);
    LangevinVelocityVerlet(Space& spc, Energy::EnergyTerm& energy, double time_step,
                           double friction_coefficient);
    LangevinVelocityVerlet(Space& spc, Energy::EnergyTerm& energy, const json& j);
    void step(PointVector& velocities, PointVector& forces) override;
    void from_json(const json& j) override;
    void to_json(json& j) const override;
};

void from_json(const json& j, IntegratorBase& i);
void to_json(json& j, const IntegratorBase& i);

class ForceMove : public Move
{
  protected:
    std::shared_ptr<IntegratorBase> integrator;
    unsigned int
        number_of_steps;    
    PointVector velocities; 
    PointVector forces;     

    size_t resizeForcesAndVelocities(); 
    void generateVelocities();          
    void _to_json(json& j) const override;
    void _from_json(const json& j) override;
    void _move(Change& change) override;
    ForceMove(Space& spc, const std::string& name, const std::string& cite,
              std::shared_ptr<IntegratorBase> integrator, unsigned int nsteps);
    ~ForceMove() override = default;

  public:
    double bias(Change&, double, double) override;
    [[nodiscard]] const PointVector& getForces() const;
    [[nodiscard]] const PointVector& getVelocities() const;
};

class LangevinDynamics : public ForceMove
{
  protected:
    LangevinDynamics(Space& spc, const std::string& name, const std::string& cite,
                     std::shared_ptr<IntegratorBase> integrator, unsigned int nsteps);

  public:
    LangevinDynamics(Space& spc, std::shared_ptr<IntegratorBase> integrator, unsigned int nsteps);
    LangevinDynamics(Space& spc, Energy::EnergyTerm&, const json&);
    LangevinDynamics(Space& spc, Energy::EnergyTerm&);
    void _to_json(json& j) const override;
    void _from_json(const json& j) override;
};

} // namespace Faunus::move