externalpotential.h

#pragma once

#include "geometry.h"
#include "group.h"
#include "aux/timers.h"
#include "aux/equidistant_table.h"
#include <set>

template <std::floating_point T> class ExprFunction;

namespace Faunus {

struct Change;
class Space;

namespace Energy {

class EnergyTerm
{
  public:
    enum class MonteCarloState
    {
        ACCEPTED,
        TRIAL,
        NONE
    };
    MonteCarloState state = MonteCarloState::NONE;
    std::string name;                                     
    std::string citation_information;                     
    TimeRelativeOfTotal<std::chrono::microseconds> timer; 
    virtual double energy(const Change& change) = 0;      
    virtual void to_json(json& j) const;                  
    virtual void sync(EnergyTerm* other_energy,
                      const Change& change); 
    virtual void init();                     
    virtual void
    updateState(const Change& change); 
    virtual void force(PointVector& forces); 
    inline virtual ~EnergyTerm() = default;
};

void to_json(json& j, const EnergyTerm& base); 

class ExternalPotential : public EnergyTerm
{
  private:
    bool act_on_mass_center = false;              
    std::set<int> molecule_ids;                   
    std::vector<std::string> molecule_names;      
    double groupEnergy(const Group& group) const; 
  protected:
    const Space& space;                                           
    std::function<double(const Particle&)> externalPotentialFunc; 
  public:
    ExternalPotential(const json& j, const Space& spc);
    double energy(const Change& j) override;
    void to_json(json& j) const override;
};

std::function<double(const Particle&)> createGouyChapmanPotential(const json& j,
                                                                  const Geometry::Chameleon&);

class CustomExternal : public ExternalPotential
{
  private:
    std::unique_ptr<ExprFunction<double>> expr;

    struct ParticleData
    { // storage for particle properties
        double charge = 0, x = 0, y = 0, z = 0;
    };

    ParticleData particle_data;
    json json_input_backup; // initial json input

  public:
    CustomExternal(const json&, Space&);
    void to_json(json&) const override;
};

class ExternalAkesson : public ExternalPotential
{
  private:
    std::string filename;                      
    bool fixed_potential = false;              
    unsigned int nstep = 0;                    
    unsigned int phi_update_interval = 0;      
    unsigned int num_rho_updates = 0;          
    unsigned int num_density_updates = 0;      
    double dielectric_constant;                
    double dz;                                 
    double bjerrum_length;                     
    double half_box_length_z;                  
    Equidistant2DTable<double> charge_profile; 
    Equidistant2DTable<double, Average<double>> rho; 
    Equidistant2DTable<double> phi;                  

    static double evalPotential(double z, double a); 
    void updateChargeDensity();                      
    void updatePotential();                          
    void saveChargeDensity();                        
    void loadChargeDensity();                        
    void to_json(json&) const override;
    void sync(EnergyTerm*, const Change&) override;

  public:
    ExternalAkesson(const json& j, const Space& spc);
    double energy(const Change& change) override;
    ~ExternalAkesson() override;
};

class Confine : public ExternalPotential
{
  public:
    enum Variant
    {
        sphere,
        cylinder,
        cuboid,
        none
    };

    Variant type = none;

  private:
    Point origo = {0.0, 0.0, 0.0};
    Point dir = {1.0, 1.0, 1.0};
    Point low;
    Point high;
    double radius = 0.0;
    double spring_constant = 0.0;
    bool scale = false;
    std::map<std::string, Variant> m = {
        {"sphere", sphere}, {"cylinder", cylinder}, {"cuboid", cuboid}};

  public:
    Confine(const json& j, Space& spc);
    void to_json(json& j) const override;
}; 

class ParticleSelfEnergy : public ExternalPotential
{
  public:
    ParticleSelfEnergy(Space&, std::function<double(const Particle&)>);
};

} // namespace Energy
} // namespace Faunus