#ifndef BERENDSEN_BAROSTAT_H
#define BERENDSEN_BAROSTAT_H

#ifndef BAROSTAT_H
#include "Barostat.h"
#endif

#ifndef MD_H
#include "MD.h"
#endif

#include "Point_3D_impl.h"
#include <math.h>

namespace MM
{

class Berendsen_barostat : public Barostat
{
    double          time_constant_;

public:
    Berendsen_barostat (Model & model, MD & simulation)
    :   Barostat (model, simulation)
    {
        Duration_unit unit;
        unit.set_ps (1.);
        time_constant_ = System_of_Units::singleton().time (unit);
    }

    void set_time_constant (Duration_unit const& Tau)
    {
        time_constant_ = System_of_Units::singleton().time (Tau);
    }
    void get_time_constant (Duration_unit & Tau)
    {
        System_of_Units::singleton().set_time (Tau, time_constant_, 1);
    }

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

    virtual void    handle_execute ()
    {
        if (model_.kit().boundary_conditions ().type_name() != "Box")
            return;

        MD_simulation_.algorithm().description ("Berendsen barostat");

        Box & box = model_.kit().boundary_conditions().box();
        //double L = box.a_mod_rectangle();

        // time step of integration algorithm
        double dt  = MD_simulation_.step_length ();

        double chi =
            1 - (dt/time_constant_) * (desired_pressure() - calc_pressure());
        if (chi <= .1)
        {
            Text message ("Berendsen barostat unstable.\n");
            message += "chi = ";
            message += chi;
            to_user().warning (message);
            return;
        }

        double chi_3 = pow (chi, 1./3.);

        if (chi_3 < 0.999 || chi_3 > 1.001)
        {
            Text message ("Berendsen barostat unstable.\n");
            message += "sqrt(chi) = ";
            message += chi_3;
            to_user().warning (message);
            if (chi_3 < 0.999) chi_3 = 0.999; 
            if (chi_3 > 1.001) chi_3 = 1.001;
        }

        box.scale_a (chi_3);
        box.scale_b (chi_3);
        box.scale_c (chi_3);

        for (int i=0;  i< model_.atom_count();  ++i)
        {
            Atom & atom = model_.atom (i);
            Point_3D_impl position (atom.position());
            position.set_x (chi_3 * position.x());
            position.set_y (chi_3 * position.y());
            position.set_z (chi_3 * position.z());
            atom.set_position (position);
        }
    }
};

}//MM

#endif //BERENDSEN_BAROSTAT_H
