#include "MD_simulator_impl.h"

#include <sstream>
//#include <string>
#include <fstream>

#include "Create.h"
#include "Program_title.h"

#include "Defs.h"
#include "User.h"
#include "System_of_Units.h"

#include "Path.h"
#include "Xyz_file.h"
#include "Own.h"
#include "Command.h"
#include "Model.h"
#include "Model_kit.h"
#include "Model_proxy.h"
#include "Project.h"
#include "Project_kit.h"
#include "Atom_group.h"
#include "Interaction.h"
#include "MD.h"
#include "MD_subject.h"
#include "Monitor.h"
#include "Atom_group_as_md_subject.h"
#include "Mass.h"
#include "Arbitrary_atom_group.h"
#include "No_interaction.h"
#include "Collect_temperature.h"
#include "Collect_molecular_mass_center_temperature.h"
#include "Collect_energy.h"
#include "Collect_potential.h"
#include "Collect_intermolecular_energy.h"
#include "Collect_density.h"
#include "Collect_density.h"
#include "Collect_pressure.h"
#include "Collect_phase_trajectory.h"

#include "PEFRL.h"
#include "VEFRL.h"
#include "OPVL.h"
#include "OVVL.h"
#include "Position_leapfrog.h"
#include "Velocity_leapfrog.h"
#include "Velocity_verlet.h"
#include "Position_verlet.h"
#include "Simulator_factory.h"
//#include "Verlet.h"

#include "Collect_phase_trajectory.h"

#include "Berendsen_thermostat.h"
#include "Andersen_thermostat.h"
#include "Lowe_thermostat.h"
#include "Nose_Hoover_thermostat.h"
#include "Isokinetic_termostat.h"

#include "Berendsen_barostat.h"

namespace MM
{
Array_of_own <MD_method>    MD_simulator_impl::method_prototype_;
MD_simulator_impl MD_simulator_impl::prototype_ (it_is_a_prototype);

MD_simulator_impl::MD_simulator_impl (Prototype &)
{
    Simulator_factory::accept ("MD", this);
    //MD_simulator::adopt (dynamic_cast <MD_simulator*> (this->clone()));
}

//MD_simulator_impl::
//MD_simulator_impl ()
//:
//    //Simulator_frame (Project::singleton().current_model()),
//
//#ifdef ASCALAPH
//    //method_     (MD_method::create ("HMC")),
//    method_     (MD_method::create ("rRESPA_4")),
//    //method_     (MD_method::create ("MTS_MD")),
//    //method_     (MD_method::create ("Shuffled")),
//    //method_     (MD_method::create ("MTS_VV")),
//    //method_     (MD_method::create ("OPVL")),
//#else
//    //method_     (MD_method::create ("MTS_VV")),
//    method_     (MD_method::create ("rRESPA_4")),
//#endif
//    //method_     (MD_method::create ("Position_verlet")),
//    //method_     (MD_method::create ("Velocity_verlet")),
//    //method_     (MD_method::create ("Velocity_leapfrog")),
//    //method_     (MD_method::create ("OPVL")),
//    //method_     (MD_method::create ("OVVL")),
//
//    //method_     (MD_method::create ("PEFRL")),
//
//    //method_     (MD_method::create ("VEFRL")),
//    simulation_ (new MD),
//    starting_pressure_   (Constant::normal_pressure),
//    final_pressure_      (Constant::normal_pressure)
//{
//    simulation_().set_method (method_()); //fix move to run
//
//    Duration_unit duration;
//#ifdef ASCALAPH
//    duration.set_step_length_ps (0.005);
//    //duration.set_step_length_ps (0.001);
//#else
//    duration.set_step_length_ps (0.005);
//#endif
//    duration.set_ps             (1000.);
//    set_duration (duration);
//
//    //set_thermostat ("Lowe");
//    set_thermostat ("Berendsen");
//    //set_thermostat ("Isokinetic");
//    //set_thermostat ("Andersen");
//    Temperature_unit t;
//    //t.set_K (50);
//    thermostat_().set_desired_temperature (t);
//
//    set_barostat ("Berendsen");
//
//
//    Model & model = Project::singleton().kit().protocol().model(); //fix to Project::singleton().curent_model()
//
//    
//    Collect_temperature * temperature_monitor = new Collect_temperature (model, "Simulator");
//    after_each_step (temperature_monitor);
//    monitor_.push_back (*temperature_monitor);
//
//    Collect_molecular_mass_center_temperature * temperature_cm_monitor = 
//        new Collect_molecular_mass_center_temperature (model, "Simulator");
//    after_each_step (temperature_cm_monitor);
//    monitor_.push_back (*temperature_cm_monitor);
//
//    Collect_density * density_monitor = new Collect_density (model, "Simulator");
//    after_each_step (density_monitor);
//    monitor_.push_back (*density_monitor);
//
//    Collect_pressure * pressure_monitor = new Collect_pressure (model, "Simulator");
//    after_each_step (pressure_monitor);
//    monitor_.push_back (*pressure_monitor);
//
//    Collect_potential * potential_monitor = new Collect_potential (model, "Simulator");
//    after_each_step (potential_monitor);
//    monitor_.push_back (*potential_monitor);
//
//    Collect_intermolecular_energy * intermolecular_monitor =
//        new  Collect_intermolecular_energy (model, "Simulator");
//    after_each_step (intermolecular_monitor);
//    monitor_.push_back (*intermolecular_monitor);
//
//    Collect_energy * energy_monitor = new Collect_energy (model, "Simulator");
//    after_each_step (energy_monitor);
//    monitor_.push_back (*energy_monitor);
//
//    Collect_phase_trajectory * phase_trajectory_monitor = 
//        new Collect_phase_trajectory (model, 
//            Path::report() +"Phase_trajectory.xmol", 6, "Simulator");//fix
//    ////////////////////////////////////////////////
//    //phase_trajectory_monitor->set_observation_period (1000/4);
//    //set_duration (Duration_unit(.1).set_step_length_ps (0.004));
//    //energy_monitor->set_observation_period (1);
//
//    //phase_trajectory_monitor->set_observation_period (1000);
//    //set_duration (Duration_unit(1).set_step_length_ps (0.001));
//    //energy_monitor->set_observation_period (10);
//
//    //// 216 H2O
//    //phase_trajectory_monitor->set_observation_period (1000/4);
//    //set_duration (Duration_unit(1).set_step_length_ps (0.004));
//    //energy_monitor->set_observation_period (10);
//
//    // 729 H2O
//    //phase_trajectory_monitor->set_observation_period (1000/4);
//    //set_duration (Duration_unit(.1).set_step_length_ps (0.0025  ));
//    //energy_monitor->set_observation_period (1);
//
//    ////////////////////////////////////////////////
//
//    after_each_step (phase_trajectory_monitor);
//    monitor_.push_back (*phase_trajectory_monitor);
//    //phase_trajectory_monitor->set_observation_period (1000);
////*/
//
//    class Save_result : public Command
//    {
//        Model & model_;
//
//    public:
//        Save_result (Model & model) : model_(model) {}
//
//        void execute ()
//        {
//            using namespace std;
//            for (int i=0;  i<1000;  ++i)
//            {
//                Text filename (Path::report());
//                filename += "step ";
//                filename += i;
//                filename += ".xmol";
//                ifstream exist (filename.c_str());
//
//                if (!exist)
//                {
//                    Xyz_file file (model_);         //fix
//                    Text old_name = file.name();
//                    file.set_name (filename);
//                    file.save (false);          
//                    file.set_name (old_name);
//
//                    //model_.save_as (filename.c_str());
//
//                    //to_user().status ();
//                    break;
//                }
//            }
//        }
//    };
//    Save_result * save_result = new Save_result (model);
//    after_stop (save_result);
//}//*/
//
MD_simulator_impl & MD_simulator_impl::
singleton ()
{
    static MD_simulator_impl instance (Model::none());
    return instance;
}

MD_simulator_impl::
MD_simulator_impl (Model & model)
:
    MD_simulator (model),

    method_     (MD_method::create ("rRESPA_4")),
    simulation_ (new MD),
    starting_pressure_   (Constant::normal_pressure),
    final_pressure_      (Constant::normal_pressure)
{
    simulation_().set_method (method_()); //fix move to run

    // fix
    // in constructor, Mass (Atom_group const & group)
    // group still undefined

    //subject_.adopt (model.kit().create_md_interface ());
    //simulation_().set_subject (subject_());

    Duration_unit duration;
#ifdef ASCALAPH
    duration.set_step_length_ps (0.005);
    //duration.set_step_length_ps (0.001);
#else
    duration.set_step_length_ps (0.005);
#endif
    duration.set_ps             (1000.);
    set_duration (duration);

    //set_thermostat ("Lowe");
    set_thermostat ("Berendsen");
    //set_thermostat ("Isokinetic");
    //set_thermostat ("Andersen");
    Temperature_unit t;
    //t.set_K (50);
    thermostat_().set_desired_temperature (t);

    set_barostat ("Berendsen");


    //Model & model = Project::singleton().kit().protocol().model(); //fix to Project::singleton().curent_model()

    Collect_temperature * temperature_monitor = new Collect_temperature (model, "Simulator");
    after_each_step (temperature_monitor);
    monitor_.push_back (*temperature_monitor);

    Collect_molecular_mass_center_temperature * temperature_cm_monitor = 
        new Collect_molecular_mass_center_temperature (model, "Simulator");
    after_each_step (temperature_cm_monitor);
    monitor_.push_back (*temperature_cm_monitor);

    Collect_density * density_monitor = new Collect_density (model, "Simulator");
    after_each_step (density_monitor);
    monitor_.push_back (*density_monitor);

    Collect_pressure * pressure_monitor = new Collect_pressure (model, "Simulator");
    after_each_step (pressure_monitor);
    monitor_.push_back (*pressure_monitor);

    Collect_potential * potential_monitor = new Collect_potential (model, "Simulator");
    after_each_step (potential_monitor);
    monitor_.push_back (*potential_monitor);

    Collect_intermolecular_energy * intermolecular_monitor =
        new  Collect_intermolecular_energy (model, "Simulator");
    after_each_step (intermolecular_monitor);
    monitor_.push_back (*intermolecular_monitor);

    Collect_energy * energy_monitor = new Collect_energy (model, "Simulator");
    after_each_step (energy_monitor);
    monitor_.push_back (*energy_monitor);

    Collect_phase_trajectory * phase_trajectory_monitor = 
        new Collect_phase_trajectory (model, 
            Path::report() +"Phase_trajectory.xmol", 6, "Simulator");//fix
    ////////////////////////////////////////////////
    //phase_trajectory_monitor->set_observation_period (1000/4);
    //set_duration (Duration_unit(.1).set_step_length_ps (0.004));
    //energy_monitor->set_observation_period (1);

    //phase_trajectory_monitor->set_observation_period (1000);
    //set_duration (Duration_unit(1).set_step_length_ps (0.001));
    //energy_monitor->set_observation_period (10);

    //// 216 H2O
    //phase_trajectory_monitor->set_observation_period (1000/4);
    //set_duration (Duration_unit(1).set_step_length_ps (0.004));
    //energy_monitor->set_observation_period (10);

    // 729 H2O
    //phase_trajectory_monitor->set_observation_period (1000/4);
    //set_duration (Duration_unit(.1).set_step_length_ps (0.0025  ));
    //energy_monitor->set_observation_period (1);

    ////////////////////////////////////////////////

    after_each_step (phase_trajectory_monitor);
    monitor_.push_back (*phase_trajectory_monitor);
    //phase_trajectory_monitor->set_observation_period (1000);
//*/

    class Save_result : public Command
    {
        Model & model_;

    public:
        Save_result (Model & model) : model_(model) {}

        void execute ()
        {
            using namespace std;
            for (int i=0;  i<1000;  ++i)
            {
                Text filename (Path::report());
                filename += "step ";
                filename += i;
                filename += ".xmol";
                ifstream exist (filename.c_str());

                if (!exist)
                {
                    Xyz_file file (model_);         //fix
                    Text old_name = file.name();
                    file.set_name (filename);
                    file.save (false);          
                    file.set_name (old_name);

                    //model_.save_as (filename.c_str());

                    //to_user().status ();
                    break;
                }
            }
        }
    };
    Save_result * save_result = new Save_result (model);
    after_stop (save_result);
}

//MD_simulator_impl::
//MD_simulator_impl (Model &      model,
//              char const * method_name)
//:
//    MD_simulator (model),
//    method_  (MD_method::create (method_name)),
//    simulation_ (new MD),
//    starting_pressure_   (Constant::normal_pressure),
//    final_pressure_      (Constant::normal_pressure)
//{
//    simulation_().set_method (method_()); //fix move to run
//}

void MD_simulator_impl::
add_comment (Text const & comment)
{
    using namespace std;
    stringstream str (comment.c_str());
    string word;

    while (str >> word)
    {
        if (word == "Berendsen" || word == "Andersen" || word == "Isokinetic")
            set_thermostat (word.c_str());
    }
}


char const * MD_simulator_impl::
method_name () const
{
    return method_().alias ();
}

void MD_simulator_impl::
set_method (char const * name)
{
    method_.adopt (MD_method::create (name));
    simulation_().set_method (method_()); //fix move to run
}


void MD_simulator_impl::
set_thermostat  (char const * name)
{
    Text type (name);

    Temperature_unit T0; 
    bool             on; 

    if (thermostat_.was_initialized ())
    {
        T0 = thermostat_().desired_temperature ();
        on = thermostat_().is_on ();
    }

    if (type == "Berendsen")
    {
        thermostat_.adopt (new Berendsen_thermostat (model(), model(), simulation_()));
    }
    else if (type == "Andersen")
    {
        thermostat_.adopt (new Andersen_thermostat (model(), model(), simulation_()));
    }
    else if (type == "Lowe")
    {
        thermostat_.adopt (new Lowe_thermostat (model(), model(), simulation_()));
    }
    else if (type == "Isokinetic")
    {
        thermostat_.adopt (new Isokinetic_termostat (model(), model(), simulation_()));
    }
    else
    {
        Text message = type;
        message += " - unknown thermostat type";
        to_user().error (message);
        return;
    }

    if (thermostat_.was_initialized ())
    {
        thermostat_().set_desired_temperature (T0);
        thermostat_().turn (on);
    }

    simulation_().set_thermostat (thermostat_());
}

void MD_simulator_impl::
set_barostat (char const * name)
{
    Text type (name);

    if (type == "Berendsen")
    {
        barostat_.adopt (new Berendsen_barostat (model(), simulation_()));
    }
    else
    {
        Text message = type;
        message += " - unknown barostat type";
        to_user().error (message);
        return;
    }

    simulation_().set_barostat (barostat_());
}

void MD_simulator_impl::
add_monitor (char const * name)
{
    Monitor * new_monitor = Monitor::create (model(), name, "Simulator");
    simulation_(). after_each_step (new_monitor);
    monitor_.push_back (*new_monitor);
}

void MD_simulator_impl::
tune ()
{
    subject_.adopt (model().kit().create_md_interface ());
    simulation_().set_subject (subject_());
    //if (thermostat_.is_valid())
    //    thermostat_().set_desired_temperature (starting_temperature_);
}

void MD_simulator_impl::
handle_start ()
{
    clear_monitors ();
    MD & md = simulation_();
    md.run();
}

int MD_simulator_impl::
steps () const
{
    return simulation_().steps ();
}

void MD_simulator_impl::
set_steps (int steps)
{
    simulation_(). set_steps (steps);
}

Duration_unit MD_simulator_impl::   // picoseconds
duration () const
{
//    return simulation_(). step_length() * Constant::time_unit_to_psec;
    Duration_unit result;
    double step_length = simulation_().step_length();
    double steps       = simulation_().steps ();
    double duration    = step_length * steps;
    System_of_Units & units = System_of_Units::singleton();
    units.set_time (result, duration, steps);

    //result.set_ps (simulation_().step_length() * Constant::time_unit_to_psec);

    //System_of_Units::singleton().set_time 
    //    (result, simulation_().step_length(), simulation_().steps ());

    return result;
}

void MD_simulator_impl::    // picoseconds
set_duration (Duration_unit duration)
{
    MD & md = simulation_();

    double step_length;
    double steps = duration.steps();

    System_of_Units & units = System_of_Units::singleton();
    units.time (duration, &step_length);
    //double  units2= new_length.ps() * Constant::psec_to_time_unit;

    md. set_step_length (step_length);
    md. set_steps ((int)(steps + .999));
}

void MD_simulator_impl::
after_each_step (Command * to_do)
{
    simulation_(). after_each_step (to_do);
}

void MD_simulator_impl::
after_stop (Command * to_do)
{
    simulation_(). after_stop (to_do);
}

void MD_simulator_impl::
eliminate_after_each_step (Command & to_do)
{
    simulation_(). eliminate_after_each_step (to_do);
}

void MD_simulator_impl::
eliminate_after_stop (Command & to_do)
{
    simulation_(). eliminate_after_stop (to_do);
}
/*
void MD_simulator::
unsubscribe (MD_monitor & to_do)
{
    algorithm(). unsubscribe (to_do);
    monitor_.    remove      (to_do);
}
//*/

void MD_simulator_impl::
update ()
{
    //fix show panel
}

}//MM
