#ifndef HMC_H
#define HMC_H

#ifndef RRESPA_4_H
#include "rRESPA_4.h"
#endif

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

#include "Thermostat.h"
#include "System_of_Units.h"

#include "Log.h"
#include "User.h"

#include <math.h>

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

//#define min(a,b) (((a) < (b)) ? (a) : (b))

namespace MM
{

class HMC : public rRESPA_4
{
protected:
    MTRand          generator_;
    int             MD_steps_;
    double          H_;
    //bool            adapt_;
    //double          adapt_up_;
    //double          adapt_down_;
    int             accept_count_, reject_count_;

public:
    explicit HMC (Prototype const &) {adopt_prototype (this);}
    explicit HMC () 
    : 
        generator_(654784),
        MD_steps_(100), 
        H_(0.),
        accept_count_(0), reject_count_(0)
        //adapt_(true), adapt_up_(1.05), adapt_down_(0.5)
    {}

    int             MD_steps () const {return MD_steps_;}
    void            set_MD_steps (double v) {MD_steps_ = (int)v;}
    
    double          accepted () const 
    {
        return double(accept_count_) / (accept_count_ + reject_count_ + 0.000001);
    }

protected:
    char const*     class_name () const {return "HMC";}
    char const*     title      () const {return "Hybrid Monte Carlo";}
    char const*     alias      () const {return "HMC";}

    MD_method *     clone () const      {return new HMC;}

    virtual void start ()
    {
        rRESPA_4::start ();
        accept_count_ = reject_count_ = 0;
    }

    virtual void step ()
    {
        setup_velocities ();
        move ();
        accept ();

        //list (MD_subject::Long2,    R2_, skin2_); 
        //F (MD_subject::Whole);
    }

    virtual void setup_velocities ()
    {
        try
        {
            simulation ().setup_velocities ();
        }
        catch(...) {to_user().fatal_error ("HMC::setup_velocities"); throw;}
    }

    virtual void move ()
    {
        try
        {
            keep_configuration ();
            H_ = H();

            for (int i=0;  i<MD_steps_;  ++i)
                rRESPA_4::step ();
        }
        catch(...) {to_user().fatal_error ("HMC::move"); throw;}
    }

    virtual void accept ()
    {
        try
        {
            Temperature_unit Tu = 
                simulation().thermostat().desired_temperature ();

            double T = System_of_Units::singleton().temperature (Tu);
            double k = System_of_Units::singleton().gas_constant();
            double beta = 1. / (k * T);

            double beta_dH = beta * (H() - H_);
            if (exp (-beta_dH) < generator_.rand())
            {
                restore_configuration ();
                //log() << "-";
                ++reject_count_;
            }
            else
            {
                //log() << "+";
                ++accept_count_;
            }
        }
        catch(...) {to_user().fatal_error ("HMC::accept"); throw;}
    }
};

}//MM

#endif //HMC_H
