#ifndef ATHMC_H
#define ATHMC_H

#ifndef HMC_H
#include "HMC.h"
#endif

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

#include "Thermostat.h"
#include "System_of_Units.h"
#include "Mach_eps.h"
#include "Energy_unit.h"

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

//#include <math.h>
//#include "MersenneTwister.h"

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

namespace MM
{

// Fischer, Cordes ans Schutte,  1998,  J Comp Chem v.19 p.1689

class ATHMC : public HMC
{
    double          K_, U_, beta_, mu_;

    double          T_plus_;
    double          c_;

    mutable int     U_count_;
    mutable double  U_mean_;
    mutable bool    auto_c_;

    mutable bool    T3_;

public:
    explicit ATHMC (Prototype const &) {adopt_prototype (this);}
    explicit ATHMC ()
    : 
        K_(0.), U_(0.), beta_(0.), mu_(0.),
        T_plus_(0.), c_ (0.), U_count_(0), U_mean_(0.), auto_c_(true)//, T3_(true)
    {
        //Temperature_unit & t = simulation().thermostat().desired_temperature();
        Temperature_unit t;     
        t.set_K (300.);
        Temperature_unit plus;  
        plus.set_K (100.);
        T_plus_ = System_of_Units::singleton().temperature (t) + 
                  System_of_Units::singleton().temperature (plus);
    }

    virtual double  U () const 
    { 
        double result = MD_method::U();
        U_mean_ += result;
        ++U_count_;
        return result;
    }

    double          T_plus () const         {return T_plus_;} 
    void            set_T_plus (double v)   {T_plus_ = v;}

    double          c () const              
    {
        //double shift = System_of_Units::singleton().energy (Energy_unit().set_kcalImol(1));
        return auto_c_ && U_count_>1 ? (U_mean_/U_count_ /*+ shift*/) : c_;} 
    void            set_c (double v)        {c_ = v;}
    
    bool            is_c_auto  () {return auto_c_;}
    void            set_c_auto (bool on) {auto_c_ = on;}

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

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

    virtual void start ()
    {
        rRESPA_4::start ();

        K_    = K();        
        U_    = U();
        beta_ = beta();     
        mu_   = mu();

        c_ = U_;
        U_mean_ = 0.;
        U_count_ = 0;

        accept_count_ = reject_count_ = 0;
    }

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

    void setup_velocities ()    
    {
        try
        {
            double k = System_of_Units::singleton().gas_constant();
            double T = 1. / (k * beta_);
            Temperature_unit Tu;
            System_of_Units::singleton().set_temperature(Tu, T);

            simulation ().setup_velocities (Tu);
        }
        catch(...) {to_user().fatal_error ("ATHMC::setup_velocities"); throw;}
    }

    void move ()
    {
        try
        {
            keep_configuration ();
            K_    = K();        
            U_    = U();
            beta_ = beta();     
            mu_   = mu();

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

    void accept ()
    {
        try
        {
            double K_new    = K();
            U_ = U();
            double beta_new = beta();
            double mu_new   = mu();
            double s        = degrees_of_freedom();

            double P_acc = (mu_new * exp(-beta_new * K_new)) / 
                           (mu_    * exp(-beta_    * K_));
            P_acc *= pow (beta_new / beta_, s/2);

            if (P_acc < generator_.rand())
            {
                restore_configuration ();
                //log() << "-";
                ++reject_count_;
            }
            else
            {
                //log() << "+";
                ++accept_count_;
            }
        }
        catch(...) {to_user().fatal_error ("ATHMC::accept"); throw;}
    }

    double beta () const
    {
        return fabs(U_ - c()) > d_mach_eps ? 
            -::log (mu()) / (U_ - c()) : 
            beta_;
    }

    double mu () const
    {
        //if (T3_)
        //    return mu_3 ();

        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_minus = 1. / (k * T);
        double beta_plus  = 1. / (k * T_plus_);
        double U_c = U_ - c();

        double result = 0.5 * (exp(-beta_minus * U_c) + exp(-beta_plus * U_c));
        return result;
    }

/*    double mu_3 () const
    {
        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_minus = 1. / (k * T);
        double beta_mid   = 1. / (k * ((T+T_plus_)/2.));
        double beta_plus  = 1. / (k * T_plus_);
        double U_c = U_ - c();

        double result = (exp(-beta_minus * U_c)
                       + exp(-beta_mid   * U_c)
                       + exp(-beta_plus  * U_c)) / 3.;
        return result;
    }//*/
};

}//MM

#endif //ATHMC_H
