#include "Test.h"

#include "Log.h"
#include "Defs.h"
#include "Path.h"

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

#include "Fine_type.h"
#include "Atom_impl.h"
#include "Element.h"
#include "Arbitrary_atom_group.h"
#include "Harmonic_potential.h"
#include "No_interaction.h"
//#include "Atom_group_configuration.h"
//#include "Velocities.h"
//#include "Force.h"
#include "Mass.h"
#include "Atom_group_as_md_subject.h"
#include "Model_impl.h"
#include "Model_kit.h"
#include "Project.h"
#include "Create.h"
#include "MD_simulator_impl.h"

#include "Collect_scalar.h"
#include "Collect_energy.h"
#include "Collect_distance.h"
#include "Collect_phase_trajectory.h"
#include "Collect_scalar_to_file.h"

//Ðóññê

namespace MM
{

class MD_test : public Test
{
    MD *link;

public:
    explicit MD_test (const Text& suite_name) : Test (suite_name) { }

    void run()
    {
        if (!Fine_type::exist ("Test type"))
            Fine_type::add ("Test type");

        //demo_consecutive_main   ();
        demo_with_MD_simulation   ();
        demo_with_model           ();                    
        demo_without_model        ();                    
        touch_algorithms          ();
        check_harmonic_oscilator  ();     
        check_harmonic_oscilator_2();     
        test_harmonic_oscilator   ();   count_tests();
        //test_bond                 ();   count_tests();

        //Fine_type::erase ("Test type");
    }


    void demo_consecutive_main   ()
    {
        Model *model = new Model_impl;
        Project::singleton().adopt (model);
        
        //======================================================

        MD_simulator_impl md(*model);
        //md.set_model (*model);
        Text file = Path::tests() + "Harmonic_H2.xyz";
        md.model().load (file);
        //model->load (file);

        Harmonic_potential *interaction 
            = new Harmonic_potential (md.model().atom(0), md.model().atom(1));
        interaction->set_equilibrium_distance (6.);
        interaction->set_force_constant       (.25);
        md.adopt_interaction (interaction);

        int steps = 20;
        md.set_duration (Duration_unit((1./steps) * 2. * __Pi__));
        md.set_steps (steps);
        //md.after_each_step ("Collect_phase_trajectory"); 
        md.add_monitor ("Collect_phase_trajectory"); 
        md.run ();

        //======================================================

         delete Project::singleton().orphan (*model);
    }

    void demo_with_MD_simulation ()
    {
        Model_impl model;
        //model.add_atom (new Atom_impl (Element ("H")));
        //model.add_atom (new Atom_impl (Element ("H")));
        Atom * atom_1 = new Atom_impl (Element ("Unknown"));
        Atom * atom_2 = new Atom_impl (Element ("Unknown"));
        atom_1->kit().set_fine_type ("Test type");
        atom_2->kit().set_fine_type ("Test type");
        model.add_atom (atom_1);
        model.add_atom (atom_2);
        model.atom(0).set_x (-4.);
        model.atom(0).set_y ( 0.);
        model.atom(0).set_z ( 0.);
        model.atom(1).set_x ( 4.);
        model.atom(1).set_y ( 0.);
        model.atom(1).set_z ( 0.);
        model.atom(0).kit().set_fine_type ("Test type");
        model.atom(1).kit().set_fine_type ("Test type");
        model.kit().interaction().set_common_on (false);
     
        Harmonic_potential *interaction 
            = new Harmonic_potential (model.atom(0), model.atom(1));
        interaction->set_equilibrium_distance (6.);
        interaction->set_force_constant       (.25);

        //MD_simulator md (model, interaction, "Position_verlet");
        model.kit().clear_interactions ();
        model.kit().adopt_interaction (interaction);
        MD_simulator_impl md (model/*, "Position_verlet"*/);
        md.set_method ("Position_verlet");
        //MD_simulator_impl md (model, "Velocity_verlet");

        Duration_unit duration;
        //duration.set_ps ((1./40.) * 2. * __Pi__ * Constant::time_unit_to_psec);
        duration.set_step_length_ps ((1./40.) * 2. * __Pi__ * Constant::time_unit_to_psec);
        duration.set_steps (40);
        md.set_duration (duration);

        //md.set_steps (40);

        //md.after_each_step (new Collect_phase_trajectory 
        //    (model, Path::bin() + "_demo_with_MD_simulation.txt", 8));
        //md.after_each_step ("Collect_phase_trajectory"); 
        //md.after_each_step ("Collect_energy");
        md.add_monitor ("Collect_energy");
        Collect_distance * distance = new Collect_distance(model, "Simulator");
        distance->set_observation_period (1);
        md.after_each_step (distance);
        
        class Collect_x_of_second_atom : public Collect_scalar
        {
        public:
            explicit Collect_x_of_second_atom (Atom_group & atom_group, Text in)
            :   Collect_scalar (atom_group, in) {}


//            MD_monitor * clone (Atom_group &atom_group, Text const& args) const
            Monitor * clone (Atom_group &atom_group, Text const&, Text in) const
                {return new Collect_x_of_second_atom (atom_group, in);}          

            Text    title     () const {return "Collect_x_of_second_atom";}
            Text    class_name() const {return "Collect_x_of_second_atom";}
            double  extract_current_value() {return atom_group().atom(1).x();}
        };

        Collect_x_of_second_atom * collect_x 
            = new Collect_x_of_second_atom (model, "Simulator");
        
        collect_x->adopt (new Collect_scalar_to_file 
            (Path::tests_rubbish() +"harmonic_oscilator___x.xls", *collect_x));

        collect_x->set_observation_period (1);

        md. after_each_step (collect_x);

        md.run ();

        DOUBLES_EQUAL (atom_1->mass(),             1.,             0.0000001);
        DOUBLES_EQUAL (distance->last_value(),     8.,             0.0001);
        DOUBLES_EQUAL (distance->average(),        6.,             0.01);       //fix!!!
        DOUBLES_EQUAL (distance->rmsd(),           2./sqrt (2.),   0.001);      //fix!!!

        double  last    = collect_x->last_value ();
        double  average = collect_x->average    ();
        double  rmsd    = collect_x->rmsd       ();
        DOUBLES_EQUAL (last,     4.,             0.0001);
        DOUBLES_EQUAL (average,  3.,             0.01);
        DOUBLES_EQUAL (rmsd,     1./sqrt (2.),   0.001); 
    }

    void demo_with_model()
    {
        //Model *model = create <Model> ();     //fix does not work
        //Own <Model> model_janitor(model);
        //model->set_MD_simulator    No!!!
        //fix to adopt:  Position_verlet md (md_subject);

        Model_impl model;
        model.add_atom (create <Atom> ());
        model.add_atom (create <Atom> ());
        model.atom(0).set_element (Element ("Unknown"));
        model.atom(1).set_element (Element ("Unknown"));
        model.atom(0).set_x (-4.);
        model.atom(1).set_x ( 4.);
     
        Harmonic_potential *interaction 
            = new Harmonic_potential (model.atom(0), model.atom(1));
        interaction->set_equilibrium_distance (6.);
        interaction->set_force_constant       (.25);
        model.kit().adopt_interaction (interaction);

        MD_subject *md_subject = model.kit().create_md_interface ();
        Own <MD_subject> subject_janitor(md_subject);

        Position_verlet method;
        MD md (*md_subject, method);
        md.set_steps (20);
        md.set_step_length ((1./20.) * 2. * __Pi__);
        md.after_each_step (new Collect_phase_trajectory 
            (model, Path::working() + "_demo_with_model.txt", 8, "Simulator"));
        md.run ();

        //model->erase ();
    }

    void demo_without_model()
    {
        Atom_impl a_1 (Element ("Unknown")); 
        Atom_impl a_2 (Element ("Unknown")); 
    
        a_1.set_x (-4.);       a_1.set_y (0.);        a_1.set_z (0.);
        a_2.set_x ( 4.);       a_2.set_y (0.);        a_2.set_z (0.);
    
        Arbitrary_atom_group model;
        model.add (a_1);
        model.add (a_2);
     
        Harmonic_potential interaction (a_1, a_2, 6., .25);
        Mass *                      masses = new Mass (model);
        Atom_group_as_md_subject    md_subject (model, interaction, masses);
        
        //Position_verlet md (md_subject);
        Position_verlet method;
        MD md (md_subject, method);
        md.set_steps (40);
        md.set_step_length ((1./40.) * 2. * __Pi__);

        Collect_distance * distance = new Collect_distance(model, "Simulator");
        distance->set_observation_period (1);
        md.after_each_step (distance);

        md. after_each_step (new Collect_phase_trajectory 
            (model, Path::working() + "_demo_without_model.txt", 8, "Simulator"));
       
        md.run ();

        DOUBLES_EQUAL (distance->last_value(),     8.,             0.001);
        DOUBLES_EQUAL (distance->average(),        6.,             0.01);
        DOUBLES_EQUAL (distance->rmsd(),           2./sqrt (2.),   0.001); 
}

    void touch_algorithms()
    {
        /*Arbitrary_atom_group     model;
        No_interaction           interaction;
        Mass *                   masses = new Mass (model);
        Atom_group_as_md_subject md_subject (model, interaction, masses);

        Velocity_leapfrog velocity_leapfrog (md_subject);
        velocity_leapfrog.set_steps (2);
        velocity_leapfrog.run ();

        Position_leapfrog position_leapfrog (md_subject);
        position_leapfrog.set_steps (2);
        position_leapfrog.run ();

        Velocity_verlet velocity_verlet (md_subject);
        velocity_verlet.set_steps (2);
        velocity_verlet.run ();

        Position_verlet position_verlet (md_subject);
        position_verlet.set_steps (2);
        position_verlet.run ();

        PEFRL pefrl (md_subject);
        pefrl.set_steps (2);
        pefrl.run ();

        VEFRL vefrl (md_subject);
        vefrl.set_steps (2);
        vefrl.run ();

        OPVL opvl (md_subject);
        opvl.set_steps (2);
        opvl.run ();

        OVVL ovvl (md_subject);
        ovvl.set_steps (2);
        ovvl.run ();
        //*/
    }

    void check_harmonic_oscilator()
    {
        // E = p*p/2 + x*x/2
        // F(x) = -x

        double  x = 1.;
        double  v = 0.; 
        int     steps = 20; //0.999659	-0.0264552
        //int     steps = 40; //0.999979  -0.00649766
        //int     steps = 100; //0.999999	-0.00103451
        //int     steps = 1000; //0.999999999946588	-1.03355e-005

        double  h = (1. / steps) * 2 * __Pi__;

        //log() << "### check_harmonic_oscilator ###\n";
        //log() << "x\tv\n";

        for (int i=0;  i<steps;  ++i)
        {
            double xI2 = x + 0.5 * h * v;
            v += -xI2 * h;
            x = xI2 + 0.5 * h * v;
            //log() << Text(x, 15) << "\t" << v << "\n";
        }

        //log() << "###  ###\n";


    }

    void check_harmonic_oscilator_2()
    {
        // E = p*p/2 + x*x/2
        // F(x) = -x

        double  x = 4.;
        double  v = 0.; 
        //int     steps = 20; // 
        int     steps = 40;  // 3.99997902017678	-0.00649766
        //int     steps = 100; //
        //int     steps = 1000; // 

        double  h = (1. / steps) * 2 * __Pi__;

        //log() << "### check_harmonic_oscilator ###\n";
        //log() << "x\tv\n";

        for (int i=0;  i<steps;  ++i)
        {
            double xI2 = x + 0.5 * h * v;   // 4 + 0.5 * 0.15707963267949 * 0
            v += -(xI2-3) * h;              // -(4-3) * 0.15707963267949 = -0.15707963267949
            x = xI2 + 0.5 * h * v;
            //log() << Text(x, 15) << "\t" << Text(v, 15) << "\n";
        }

        //log() << "###  ###\n";


    }

    void test_harmonic_oscilator()
    {
        Atom_impl a_1 (Element ("Unknown")); 
        Atom_impl a_2 (Element ("Unknown")); 
    
        a_1.kit().set_fine_type ("Test type");
        a_2.kit().set_fine_type ("Test type");
    
        a_1.set_x (-4.);       a_1.set_y (0.);        a_1.set_z (0.);
        a_2.set_x ( 4.);       a_2.set_y (0.);        a_2.set_z (0.);
    
        Arbitrary_atom_group model;
        //Atom_groupe_for_simulation_impl model;    ?= Interaction
        //model.adopt (Interaction);
        //Atom_groupe_for_simulation::update_force()
        model.add (a_1);
        model.add (a_2);

        Harmonic_potential potential (a_1, a_2, 6., .25);

        TEST ("Zero atoms velocities.", 
                a_1.kit().velocity().length() == 0 &&
                a_2.kit().velocity().length() == 0);

        //Atom_group_configuration    coordinats (model, potential);
        //Velocities                  velocities (model);
        //Force                       forces     (model);
        Mass *                      masses = new Mass (model);

        Atom_group_as_md_subject    md_subject (model, potential, masses);
        //coordinats.after_change_update (interactions); //fix remove

        int steps = 40;

        // 2 steps
        //OVVL method;               // 4.19657831400307	0.12755976746791
        //OPVL method;               // 4.19657831400307	-3.38507720822752
        //VEFRL method;              // 3.536759157688  	0.37653802200109
        //PEFRL method;              // 4.08545871585761	-0.353089008708788

        // 3 steps
        //Velocity_verlet method;    //-0.21620039027293	-5.33300780023111
        //Position_verlet method;    // 4.86981526851118	20.1556343558176
        //OVVL method;               // 3.90176348547584	-0.493551681502117
        //OPVL method;               // 3.90176348547584	-0.378526957286145
        //VEFRL method;              // 3.98792758162415	0.153372306638364
        //PEFRL method;              // 3.99998630969858	0.00497659413378559     #

        // 4 steps
        //Velocity_verlet method;    // 3.58693563637912	-0.501155821751533
        //Position_verlet method;    // 3.58693563637912	-1.3079895120389
        //OVVL method;               // 3.977358664623  	-0.216435975791591
        //OPVL method;               // 3.977358664623  	-0.206851197092386
        //VEFRL method;              // 3.99914775322368	0.0413288778002163
        //PEFRL method;              // 3.99999774802763	0.00209733237596432     #

        // 5 steps
        //Velocity_verlet method;    // 3.87239574605102	-0.380265159535209
        //Position_verlet method;    // 3.87239574605102	-0.628313313173695
        //OVVL method;               // 3.99166037598745	-0.129759995682101
        //OPVL method;               // 3.99166037598745	-0.128003230958131
        //VEFRL method;              // 3.99987739522267	0.0156733953811272
        //PEFRL method;              // 3.99999958463507	0.000907517985479278    #

        // 10 steps
        ////Velocity_leapfrog method;  //-  
        ////Position_leapfrog method;  //-   
        //Velocity_verlet method;    // 3.99414844241952	-0.102553413839177
        //Position_verlet method;    // 3.99414844241952	-0.113783383682848
        //OVVL method;               // 3.99953816935925	-0.0303870420759327
        //OPVL method;               // 3.99953816935925	-0.0303895322106023
        //VEFRL method;              // 3.99999960549692	0.000888330725603279
        //PEFRL method;              // 3.99999999826429	5.89049682861442e-005

        // 20 steps
        //Velocity_verlet method;    // 3.99965863768855	-0.025802430409389
        //Position_verlet method;    // 3.99965863768855	-0.0264551859588144
        //OVVL method;               // 3.99997191917206	-0.00749334568150172
        //OPVL method;               // 3.99997191917206	-0.00749476532082144
        //VEFRL method;              // 3.99999999852889	5.42425332220389e-005
        //PEFRL method;              // 3.99999999999315	3.7002480742232e-006

        // 40 steps
        Velocity_verlet method;    // 3.99997902017678	-0.006457583202051
        //Position_verlet method;    // 3.99997902017678	-0.00649766406068181
        //OVVL method;               // 3.99999825673411	-0.00186717136988901
        //OPVL method;               // 3.99999825673411	-0.00186727838525214
        //VEFRL method;              // 3.99999999999432	3.37073729478679e-006
        //PEFRL method;              // 3.99999999999997	2.31507591860747e-007

        // 100 steps
        //Velocity_verlet method;    // 3.99999946542013	-0.00103349128579949
        //Position_verlet method;    // 3.99999946542013	-0.0010345123085212
        //OVVL method;               // 3.99999995545607	-0.000298474446218609
        //OPVL method;               // 3.99999995545607	-0.000298477314373276
        //VEFRL method;              // 4               	8.61526558409292e-008
        //PEFRL method;              // 4               	5.92824694889416e-009

        // 1000 steps
        //Velocity_verlet method;    // 3.99999999994658	-1.03354204591837e-005
        //Position_verlet method;    // 3.99999999994659	-1.03355224720233e-005
        //OVVL method;               // 3.99999999999555	-2.98423127052696e-006
        //OPVL method;               // 3.99999999999556	-2.98423156074598e-006
        //VEFRL method;              // 4               	8.62063740704466e-012
        //PEFRL method;              // 4.00000000000001	5.97753617437302e-013

        MD md (md_subject, method);

        //md.set_steps (10*steps);
        md.set_steps (steps);
        md.set_step_length ((1./steps) * 2. * __Pi__);
        //md.set_thermostat (...);
        //md.set_barostat   (...);
        //md.after_each_step_calculate (update_veiws);
        //md.after_each_step_calculate (configuration_trajectory);

        class Collect_x_of_second_atom : public Collect_scalar
        {
        public:
            explicit Collect_x_of_second_atom (Atom_group & atom_group, Text in)
            :   Collect_scalar (atom_group, in) {}


//            MD_monitor * clone (Atom_group &atom_group, Text const& args) const
            Monitor * clone (Atom_group &atom_group, Text const&, Text in) const
                {return new Collect_x_of_second_atom (atom_group, in);}          

            Text    title     () const {return "Collect_x_of_second_atom";}
            Text    class_name() const {return "Collect_x_of_second_atom";}
            double  extract_current_value() {return atom_group().atom(1).x();}
        };

        Collect_x_of_second_atom * collect_x 
            = new Collect_x_of_second_atom (model, "Simulator");
        
        collect_x->adopt (new Collect_scalar_to_file 
            (Path::tests_rubbish() +"harmonic_oscilator___x.xls", *collect_x));

        collect_x->set_observation_period (1);

        md. after_each_step (collect_x);
        


        md. after_each_step (new Collect_phase_trajectory 
            (model, Path::tests_rubbish() 
            + "harmonic_oscilator___phase_trajectory.xls", 15, "Simulator"));


        
        Collect_energy * collect_energy = new Collect_energy (model, "Simulator");
        collect_energy->adopt (new Collect_scalar_to_file (Path::tests_rubbish() 
            + "harmonic_oscilator___energy.xls", *collect_energy));
        md. after_each_step (collect_energy);



        Collect_distance * collect_distance = new Collect_distance (model, "Simulator");
        collect_distance->adopt (new Collect_scalar_to_file (Path::tests_rubbish() 
            + "harmonic_oscilator___distance.xls", *collect_distance));
        md. after_each_step (collect_distance);

        //md.after_each_step_calculate (gradient);
        md.run ();//*/

        double  last    = collect_x->last_value ();
        double  average = collect_x->average    ();
        double  rmsd    = collect_x->rmsd       ();
        DOUBLES_EQUAL (last,     4.,             0.0001);
        DOUBLES_EQUAL (average,  3.,             0.01);
        DOUBLES_EQUAL (rmsd,     1./sqrt (2.),   0.001); 
        //0.70710678118654752440084436210485
    }

    void test_bond ()
    {
        Model_impl model;
        Text file = Path::tests() + "quadratic_potential.mlm";
        model.load (file, "auto detect");

        MD_simulator_impl md (model);
        Duration_unit duration;
        duration.set_step_length_ps (0.001);
        duration.set_steps (1000);
        md.set_duration (duration);

        DOUBLES_EQUAL (model.atom(0).position().x(),  -4.000,  0.01);
        md.run ();
        DOUBLES_EQUAL (model.atom(0).position().x(),  -2.965,  0.01);
    }

};

MD_test test_md ("correctness");
//MD_test test_md ("demo");

}//MM
