#ifndef LOWE_THERMOSTAT_H
#define LOWE_THERMOSTAT_H

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

#ifndef VECTOR_3D_IMPL_H
#include "Vector_3D_impl.h"
#endif

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

#ifndef POINT_3D_H
#include "Point_3D.h"
#endif

#ifndef DBC_H
#include "DbC.h"
#endif

#ifdef WIN32
#pragma warning(disable : 4800)
#pragma warning(disable : 4996)
#endif
#include "MersenneTwister.h"

#include "anint.h"

namespace MM
{

// Koopman and Lowe (2006) J. Chem. Phis. v.124. p.204103

class Lowe_thermostat : public Thermostat, protected DbC
{
    double          nu_;
    MTRand          distr_;

public:
    Lowe_thermostat (Model & model, Atom_group & atom_group, MD & simulation)
    :   
        Thermostat (model, atom_group, simulation), nu_(0.1), distr_(1966) {}

    char const* title () const      {return "Lowe";}
    double              response_constant () const       {return nu_;}
    void                set_response_constant (double v) {nu_ = v;}

    virtual void    handle_execute ()
    {
        //Hydrogen_atoms group (atom_group_);
        Atom_group & group = atom_group_;
        double LI2_2 = 1000;
        bool    box = false;
        double  L = 0., lIL = 0.;

        MD_simulation_.algorithm().description ("Lowe thermostat");

        if (model_.kit().boundary_conditions ().is_a_box())
        {
            box = true;
            L = model_.kit().boundary_conditions ().box ().a_mod_rectangle();
            lIL = 1. / L;
            LI2_2 = L*L/4.;
        }

        double dt    = MD_simulation_.step_length ();                           // time step of integration algorithm
        double T0    = T0_ * System_of_Units::singleton().gas_constant();
        int    count = group.atom_count();

        for (int i=0;  i<count;  ++i)
        {
            if (distr_.rand() < nu_*dt)
            {
                Atom & atom_i       = group.atom (i);
                Vector_3D_impl v_i  (atom_i.kit().velocity ());
                double mass_i       = atom_i.mass();
                //Point_3D const& r_i = atom_i.position();
                double x_i = atom_i.x();
                double y_i = atom_i.y();
                double z_i = atom_i.z();

                int j = distr_.randInt (count-1);
                if (i == j)
                    continue;
                CHECK ("j >= 0", j >= 0);
                CHECK ("j < 0",  j < count);

                Atom & atom_j       = group.atom (j);
                Vector_3D_impl v_j  (atom_j.kit().velocity ());
                double mass_j       = atom_j.mass();
                //Point_3D const& r_j = atom_j.position();
                double x_j = atom_j.x();
                double y_j = atom_j.y();
                double z_j = atom_j.z();
                
                double x = x_i - x_j;
                double y = y_i - y_j ;
                double z = z_i - z_j ;
                
                if (box)
                {
                    x -= L * anint (x * lIL);
                    y -= L * anint (y * lIL);
                    z -= L * anint (z * lIL);
                }

                Vector_3D_impl sigma (x,y,z);
                if (sigma.length_2() > LI2_2)
                    continue;
                sigma.norm();

                double mu = (mass_i * mass_j) / (mass_i + mass_j);
                CHECK ("mu > 0.", mu > 0.);

                double lambda = gauss() * sqrt (T0 / mu);

                Vector_3D_impl v(v_i);  v -= v_j;

                sigma *= lambda - v.dot(sigma);

                Vector_3D_impl k_i (sigma), k_j (sigma);
                k_i *= mu / mass_i;
                k_j *= mu / mass_j;
                CHECK ("k_i < 100.", k_i.length() < 10000.);
                CHECK ("k_j < 100.", k_j.length() < 10000.);

                v_i += k_i;
                v_j -= k_j;

                atom_i.kit().set_velocity (v_i);
                atom_j.kit().set_velocity (v_j);
            }
        }
    }

    double gauss ()    {return distr_.randNorm ();}
};

}//MM

#endif //LOWE_THERMOSTAT_H
