#ifndef BERENDSEN_THERMOSTAT_H
#define BERENDSEN_THERMOSTAT_H

#ifndef THERMOSTAT_H
#include "Thermostat.h"
#endif

#ifndef ATOM_KIT_H
#include "Atom_kit.h"
#endif

#ifndef HYDROGEN_ATOMS_H
#include "Hydrogen_atoms.h"
#endif

#ifndef INITIALIZATION_STAGE_H
#include "Initialization_stage.h"
#endif

#ifndef  MACH_EPS_H
#include "Mach_eps.h"
#endif

namespace MM
{

class Berendsen_thermostat : public Thermostat
{
    double          Tau_;  // time constant
    int             count_;
    int             renew_period_;

public:
    Berendsen_thermostat (Model & model, Atom_group & atom_group, MD & simulation)
    :
        Thermostat (model, atom_group, simulation)
    {
        Duration_unit tau;
        tau.set_ps (0.4);
        //tau.set_ps (0.1);
        //tau.set_ps (0.04);
        Tau_ = System_of_Units::singleton().time (tau);
        count_ = 0;
        renew_period_ = 1000;
    }

    char const* title () const              {return "Berendsen";}

    int         renew_period () const       {return renew_period_;}
    void        set_renew_period (int v)    {renew_period_ = v;}

    virtual void    handle_execute ()
    {
        MD_simulation_.algorithm().description ("Berendsen thermostat");
        //Hydrogen_atoms group (atom_group_);
        Atom_group & group = atom_group_;

        if (++count_ % renew_period_) //fix
        {
            double dt  = MD_simulation_.step_length ();  // time step of integration algorithm
            double T   = calc_temperature (group);   // actual temperature
    
            double lambda = sqrt (1 + dt/Tau_ * (T0_/(T+d_mach_eps_2) - 1));
    
    
            int    count = group.atom_count();
    
            for (int i=0;  i<count;  ++i)
            {
                Atom & atom = group.atom (i);
                atom.kit().scale_velocity (lambda);
            }
        }
        else
        {
            Temperature_unit t;
            System_of_Units::singleton().set_temperature (t, T0_);
            Initialization_stage::setup_velocities (group, t);
        }
    }

    double response_constant () const // ps
    {
        Duration_unit tau;
        System_of_Units::singleton().set_time (tau, Tau_, 1);
        return tau.ps ();
    }
    void set_response_constant (double v)  // ps
    {
        Duration_unit tau;
        tau.set_ps (v);
        Tau_ = System_of_Units::singleton().time (tau);
    }
};

}//MM

#endif //BERENDSEN_THERMOSTAT_H
