#include "Model_as_md_subject.h"

#include "Program_title.h"
#include "Model.h"
#include "Model_kit.h"
#include "Atom.h"
#include "Atom_kit.h"
#include "Atom_cache.h"
#include "Atom_group.h"
#include "Interaction.h"
#include "Mass.h"
#include "Vector_3D_impl.h"
#include "Initialization_stage.h"

#include "Log.h"

#include "Profiler.h"
#include "Timer.h"
#define TIME(ref) ;

//#include <cmnintrin.h>
//#include <xmmintrin.h>

namespace MM
{

Model_as_md_subject::
Model_as_md_subject (Model & model)
:
    model_(model),// interaction_(interaction),
    masses_(new Mass (model)),
    //atom_group_as_md_subject_(model_, interaction_, new Mass(model))
    first_bad_iteration_  (true),
#ifndef ASCALAPH
    lIforce_limit_  (1. / 1000.),
#else
    lIforce_limit_  (1. / 300.),
    //lIforce_limit_  (1. / 100),
#endif
    H_monitor_(model, "md_subject"),
    K_monitor_(model, "md_subject"),
    U_monitor_(model, "md_subject")
{
}

//void Model_as_md_subject::
//F (Hint hint)
//{
//    update_interactions (hint);
//}

void Model_as_md_subject::start  (){model_.kit().interaction().start();}
void Model_as_md_subject::finish (){model_.kit().interaction().finish();}
void Model_as_md_subject::zero_virial ()
                                   {model_.kit().interaction().zero_virial();}

void Model_as_md_subject::
list (Hint hint, double skin)
{
    try
    {
        switch (hint)
        {
        //case Long : model_.kit().interaction().list (chunk3_, skin); break;
        //case Short: model_.kit().interaction().list (chunk1_, skin); break;//fix

        case Whole: model_.kit().interaction().list (skin); break;

        case Long :
        case Long3: model_.kit().interaction().list (chunk3_, skin); break;
        //case Long2: model_.kit().interaction().list (chunk2_, skin); break;
        //case Long1: model_.kit().interaction().list (chunk1_, skin); break;
        default:
            FLAW ("Model_as_md_subject::list - Wrong hint");
        }
    }
    catch(...) 
    {
        Text message ("Model_as_md_subject::list (Hint ");
        message += hint;
        message += ", double skin=";
        message += skin;
        message += ")";
        to_user().fatal_error (message);
        throw;
    }
}

void Model_as_md_subject::
list (Hint hint, double R, double skin)
{
    try
    {
        switch (hint)
        {
        //case Long : model_.kit().interaction().list (chunk3_, R+skin);          break;
        case Short: model_.kit().interaction().list (chunk1_, R+skin, chunk3_); break;//fix

        //case Long3: model_.kit().interaction().list (chunk3_, R+skin);          break;
        case Long2: model_.kit().interaction().list (chunk2_, R+skin, chunk3_); break;
        case Long1: model_.kit().interaction().list (chunk1_, R+skin, chunk2_); break;
        default:
            FLAW ("Model_as_md_subject::list - Wrong hint");
        }
    }
    catch(...) 
    {
        Text message ("Model_as_md_subject::list (Hint ");
        message += hint;
        message += ", double R=";
        message += R;
        message += ", double skin=";
        message += skin;
        message += ")";
        to_user().fatal_error (message);
        throw;
    }
}

double Model_as_md_subject::
degrees_of_freedom ()
{
    return 3. * model_.atom_count();                //fix
}

double Model_as_md_subject::
H () 
{
    try
    {
        return H_monitor_.extract_current_value ();
    }
    catch(...) 
    {
        Text message ("Model_as_md_subject::H");
        to_user().fatal_error (message);
        throw;
    }
}

double Model_as_md_subject::
K () 
{
    try
    {
        return K_monitor_.extract_current_value ();
    }
    catch(...) 
    {
        Text message ("Model_as_md_subject::K");
        to_user().fatal_error (message);
        throw;
    }
}

double Model_as_md_subject::
U ()
{
    try
    {
        return U_monitor_.extract_current_value ();
    }
    catch(...) 
    {
        Text message ("Model_as_md_subject::U");
        to_user().fatal_error (message);
        throw;
    }
}

void Model_as_md_subject::
F (Hint hint)
{
    try
    {
        //int stubs = 3; //fix

        //for (int i=0;  i<size;  ++i)
        //    model_.atom(i).kit().set_potential_force (Vector_3D_impl (0,0,0));

        bool ad_hoc = model_.kit().interaction().common().ad_hoc();
        model_.kit().force_field().start_work (!ad_hoc);
        {
            switch (hint)
            {
            case Whole:
                model_.kit().interaction().add_force ();
                break;

            case Long:
            case Long3:
                model_.kit().interaction().add_long_force (chunk3_);
                extract_long (Long3);
                break;
            //case Long2:
            //    model_.kit().interaction().add_long_force (chunk2_);
            //    extract_long (Long2);
            //    break;
            //case Long1:
            //    model_.kit().interaction().add_long_force (chunk1_);
            //    extract_long (Long1);
            //    break;
            case Short:
                TIME ("Short <              ");
                model_.kit().interaction().add_short_force ();
                TIME ("Short >              ");
                break;
            default:
                FLAW ("Model_as_md_subject::F - Wrong hint");
            }
        }
        if  (!ad_hoc && model_.kit().force_field().was_problems ())
            throw Text ("Force Field problem");

        TIME ("cut_large_forces <              ");
        cut_large_forces ();
        TIME ("cut_large_forces >              ");
    }
    catch(...) 
    {
        Text message ("Model_as_md_subject::F (Hint ");
        message += hint;
        message += ")";
        to_user().fatal_error (message);
        throw;
    }
}

void Model_as_md_subject::
F (Hint hint, double R, double Switching)
{
    try
    {
        //int stubs = 3; //fix
        //int size = model_.atom_count() /*+ stubs*/;

        //for (int i=0;  i<size;  ++i)
        //    model_.atom(i).kit().set_potential_force (Vector_3D_impl (0,0,0));

        bool ad_hoc = model_.kit().interaction().common().ad_hoc();
        model_.kit().force_field().start_work (!ad_hoc);
        {
            switch (hint)
            {
            case Whole:
                model_.kit().interaction().add_force ();
                break;
            //case Long:
            //    TIME (" { Long ");
            //    f_long3_x_.assign(size, 0.); 
            //    f_long3_y_.assign(size, 0.); 
            //    f_long3_z_.assign(size, 0.);
            //    model_.kit().interaction().add_long_force 
            //        (chunk3_);
            //    TIME (" } Long ");
            //    break;

            //case Long3:
            //    model_.kit().interaction().add_long_force (R, Switching, chunk3_);
            //    extract_long (Long3);
            //    break;
            case Long2:
                model_.kit().interaction().add_long_force (R, Switching, chunk2_);
                extract_long (Long2);
                break;
            case Long1:
                model_.kit().interaction().add_long_force (R, Switching, chunk1_);
                extract_long (Long1);
                break;
            //case Short:
            //    model_.kit().interaction().add_short_force ();
            //    break;
            default:
                FLAW ("Model_as_md_subject::F - Wrong hint");
            }
        }
        if  (!ad_hoc && model_.kit().force_field().was_problems ())
            throw Text ("Force Field problem");

        cut_large_forces ();
    }
    catch(...) 
    {
        Text message ("F (Hint ");
        message += hint;
        message += ", double R=";
        message += R;
        message += ", double Switching=";
        message += Switching;
        message += ")";
        to_user().fatal_error (message);
        throw;
    }
}

void Model_as_md_subject::
clear_moment (Hint hint)
{
    int s = model_.atom_count();

    switch (hint)
    {
    case Long1:
        {
            sum1_x_.clear ();  
            sum1_y_.clear ();
            sum1_z_.clear ();
            sum1_x_.resize (s, 0.);  
            sum1_y_.resize (s, 0.);
            sum1_z_.resize (s, 0.);
        }
        break;
    }
}

void Model_as_md_subject::
kick (double dT, Hint hint)
{
    profile().interactions___add_force___others ();

    int i, size = model_.atom_count();

    switch (hint)
    {
    case Whole:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl (cache.force_x___, 
                                  cache.force_y___, 
                                  cache.force_z___);
            //dV += atom.kit().potential_force();

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;

    case Long:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            //Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl(f_long3_x_[i], f_long3_y_[i], f_long3_z_[i]);
            dV -= Vector_3D_impl(f_long_x_ [i], f_long_y_ [i], f_long_z_ [i]);

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;

    case Long3:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            //Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl(f_long3_x_[i], f_long3_y_[i], f_long3_z_[i]);
            dV -= Vector_3D_impl(f_long2_x_[i], f_long2_y_[i], f_long2_z_[i]);

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;
    case Long2:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            //Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl(f_long2_x_[i], f_long2_y_[i], f_long2_z_[i]);
            dV -= Vector_3D_impl(f_long_x_[i],  f_long_y_[i],  f_long_z_[i]);

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;
    case Long1:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            //Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl(f_long_x_[i],  f_long_y_[i],  f_long_z_[i]);
            //dV += atom.kit().nonbonded_force();

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;
    case Short:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
    
            //_prefetch (&model_.atom(i+1));

            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl 
                (cache.F_bond_stretch_x___ + cache.F_angle_x___ + 
                 cache.F_torsion_x___      + cache.F_improper_x___,

                 cache.F_bond_stretch_y___ + cache.F_angle_y___ + 
                 cache.F_torsion_y___      + cache.F_improper_y___, 

                 cache.F_bond_stretch_z___ + cache.F_angle_z___ + 
                 cache.F_torsion_z___      + cache.F_improper_z___);

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;
    }


    //for (int i=0;  i<size;  ++i)
    //{
    //    // v = v0 + dT * F / m
    //    Atom &       atom = model_.atom(i);
    //    Vector_3D_impl  dV (0,0,0);
    //    Atom_cache & cache = *atom.kit().cache___;

    //    switch (hint)
    //    {
    //    case Whole:
    //        dV += Vector_3D_impl (cache.force_x___, 
    //                              cache.force_y___, 
    //                              cache.force_z___);
    //        //dV += atom.kit().potential_force();
    //        break;

    //    case Long:
    //        dV += Vector_3D_impl(f_long3_x_[i], f_long3_y_[i], f_long3_z_[i]);
    //        dV -= Vector_3D_impl(f_long_x_ [i], f_long_y_ [i], f_long_z_ [i]);
    //        break;

    //    case Long3:
    //        dV += Vector_3D_impl(f_long3_x_[i], f_long3_y_[i], f_long3_z_[i]);
    //        dV -= Vector_3D_impl(f_long2_x_[i], f_long2_y_[i], f_long2_z_[i]);
    //        break;
    //    case Long2:
    //        dV += Vector_3D_impl(f_long2_x_[i], f_long2_y_[i], f_long2_z_[i]);
    //        dV -= Vector_3D_impl(f_long_x_[i],  f_long_y_[i],  f_long_z_[i]);
    //        break;
    //    case Long1:
    //        dV += Vector_3D_impl(f_long_x_[i],  f_long_y_[i],  f_long_z_[i]);
    //        //dV += atom.kit().nonbonded_force();
    //        break;
    //    case Short:
    //        dV += Vector_3D_impl (cache.F_bond_stretch_x___ +
    //                              cache.F_angle_x___ + 
    //                              cache.F_torsion_x___ +
    //                              cache.F_improper_x___,

    //                              cache.F_bond_stretch_y___ +
    //                              cache.F_angle_y___ + 
    //                              cache.F_torsion_y___ +
    //                              cache.F_improper_y___, 

    //                              cache.F_bond_stretch_z___ +
    //                              cache.F_angle_z___ + 
    //                              cache.F_torsion_z___ +
    //                              cache.F_improper_z___);
    //        break;
    //    }

    //    if (atom.kit().was_frozen())
    //        continue;

    //    dV *= dT * masses_().inverse(i);

    //    atom.kit().add_velocity (dV);
    //}

    profile().kick ();
}

void Model_as_md_subject::
kick (int substeps, double dT, Hint hint)
{
    profile().interactions___add_force___others ();

    double Ilsubsteps = 1. / substeps;
    int i, size = model_.atom_count();

    switch (hint)
    {
    case Whole:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl (cache.force_x___, 
                                  cache.force_y___, 
                                  cache.force_z___);
            //dV += atom.kit().potential_force();

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;

    case Long:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            //Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl (f_long3_x_[i], f_long3_y_[i], f_long3_z_[i]);
            //dV -= Vector_3D_impl(f_long_x_ [i], f_long_y_ [i], f_long_z_ [i]);
            dV -= Vector_3D_impl (sum1_x_[i] * Ilsubsteps, 
                                  sum1_y_[i] * Ilsubsteps, 
                                  sum1_z_[i] * Ilsubsteps);

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;

    case Long3:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            //Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl(f_long3_x_[i], f_long3_y_[i], f_long3_z_[i]);
            dV -= Vector_3D_impl(f_long2_x_[i], f_long2_y_[i], f_long2_z_[i]);

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;
    case Long2:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            //Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl(f_long2_x_[i], f_long2_y_[i], f_long2_z_[i]);
            dV -= Vector_3D_impl(f_long_x_[i],  f_long_y_[i],  f_long_z_[i]);

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;
    case Long1:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            //Atom_cache & cache = *atom.kit().cache___;

            dV += Vector_3D_impl(f_long_x_[i],  f_long_y_[i],  f_long_z_[i]);
            //dV += atom.kit().nonbonded_force();

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;
    case Short:
        for (i=0;  i<size;  ++i)
        {
            // v = v0 + dT * F / m
    
            //_mm_prefetch ((char *)&model_.atom(i+3), _MM_HINT_NTA);

            Atom &       atom = model_.atom(i);
            Vector_3D_impl  dV (0,0,0);
            Atom_cache & cache = *atom.kit().cache___;

            //_mm_prefetch ((char *)(&cache+1), _MM_HINT_NTA);

            dV += Vector_3D_impl 
                (cache.F_bond_stretch_x___ + cache.F_angle_x___ + 
                 cache.F_torsion_x___      + cache.F_improper_x___,

                 cache.F_bond_stretch_y___ + cache.F_angle_y___ + 
                 cache.F_torsion_y___      + cache.F_improper_y___, 

                 cache.F_bond_stretch_z___ + cache.F_angle_z___ + 
                 cache.F_torsion_z___      + cache.F_improper_z___);

            if (atom.kit().was_frozen())
                continue;

            dV *= dT * masses_().inverse(i);
            atom.kit().add_velocity (dV);
        }
        break;
    }

    profile().kick ();
}

void Model_as_md_subject::
extract_long (Hint hint)
{
    int i, s = model_.atom_count() /*+ stubs*/;
    
    switch (hint)
    {
        case Long3:
        case Long:
            f_long3_x_.resize(s); f_long3_y_.resize(s); f_long3_z_.resize(s);

            for (i=0;  i<s;  ++i)
            {
                Atom &       atom  = model_.atom(i);
                Atom_cache & cache = *atom.kit().cache___;
                f_long3_x_[i] = cache.F_nonbonded_x___;
                f_long3_y_[i] = cache.F_nonbonded_y___;
                f_long3_z_[i] = cache.F_nonbonded_z___;
            }
            break;
        case Long2:
            f_long2_x_.resize(s); f_long2_y_.resize(s); f_long2_z_.resize(s);

            for (i=0;  i<s;  ++i)
            {
                Atom &       atom  = model_.atom(i);
                Atom_cache & cache = *atom.kit().cache___;
                f_long2_x_[i] = cache.F_nonbonded_x___;
                f_long2_y_[i] = cache.F_nonbonded_y___;
                f_long2_z_[i] = cache.F_nonbonded_z___;
            }
            break;
        case Long1:
            f_long_x_.resize(s); f_long_y_.resize(s); f_long_z_.resize(s);

            for (i=0;  i<s;  ++i)
            {
                Atom &       atom  = model_.atom(i);
                Atom_cache & cache = *atom.kit().cache___;
                f_long_x_[i] = cache.F_nonbonded_x___;
                f_long_y_[i] = cache.F_nonbonded_y___;
                f_long_z_[i] = cache.F_nonbonded_z___;
                //sum1_x_[i] += f_long_x_[i] = cache.F_nonbonded_x___;
                //sum1_y_[i] += f_long_y_[i] = cache.F_nonbonded_y___;
                //sum1_z_[i] += f_long_z_[i] = cache.F_nonbonded_z___;
            }
            break;
        default:
            FLAW ("Model_as_md_subject::extract_long - Wrong hint");
    }
}

void Model_as_md_subject::
drift (double dT)
{
    profile().interactions___add_force___others ();

    int size = model_.atom_count();

    for (int i=0;  i<size;  ++i)
    {
        // x = x0 + dT * v0
        Atom &          atom = model_.atom(i);

        if (atom.kit().was_frozen())
            continue;

        Vector_3D_impl  movement (atom.kit().velocity());
        movement *= dT;
        //atom.move (movement, true);

        Atom_cache & cache = *atom.kit().cache___;
        cache.x___ += movement.x();
        cache.y___ += movement.y();
        cache.z___ += movement.z();
    }
    //update_interactions ();
    profile().drift ();
}

void Model_as_md_subject::
keep_configuration ()
{
    int size = model_.atom_count();
    
    configuration_.clear();
    //configuration_.reserve (size);

    for (int i=0;  i<size;  ++i)
        configuration_.push_back (model_.atom(i).position());
}

void Model_as_md_subject::
restore_configuration ()
{
    int size = model_.atom_count();

    for (int i=0;  i<size;  ++i)
        model_.atom(i).set_position (configuration_[i]);
}

void Model_as_md_subject::
update (bool force)
{
    profile().interactions___add_force___others ();

    int size = model_.atom_count();

    if (model_.kit().visual().is_on() || force)
        for (int i=0;  i<size;  ++i)
        {
            Atom & atom = model_.atom(i);
            atom.update ();
        }

    model_.kit().event().conformation_changed ();

    profile().update ();
}

bool Model_as_md_subject::locked () const
{
    return model_.kit().interaction().locked();
}
void Model_as_md_subject::set_lock (bool on)
{
    model_.kit().interaction().set_lock (on);
}


void Model_as_md_subject::
setup_velocities (Temperature_unit const& t)
{
    Initialization_stage::setup_velocities (model_, t);
}

//void Model_as_md_subject::
//update_interactions (Hint hint)
//{
//    //fix describe
//    //fix dynamic check or DbC
//    int size = model_.atom_count();
//
//    for (int i=0;  i<size;  ++i)
//        model_.atom(i).kit().set_potential_force (Vector_3D_impl (0,0,0));
//
//    bool ad_hoc = model_.kit().interaction().common().ad_hoc();
//    model_.kit().force_field().start_work (!ad_hoc);
//    {
//        switch (hint)
//        {
//        case Whole:
//            model_.kit().interaction().add_force ();
//            break;
//        //case Long:
//        //    model_.kit().interaction().add_long_force ();
//        //    break;
//        case Short:
//            model_.kit().interaction().add_short_force ();
//            break;
//        }
//    }
//    if  (!ad_hoc && model_.kit().force_field().was_problems ())
//        throw Text ("Force Field problem");
//
//    cut_large_forces ();
//}

void Model_as_md_subject::
cut_large_forces ()
{
    if (!to_cut_large_forces ())
        return;

    bool log_it = first_bad_iteration_;
    int size = model_.atom_count();
    double lIforce_limit_2 = lIforce_limit_*lIforce_limit_;

    for (int i=0;  i<size;  ++i)
    {
        Atom &         atom  = model_.atom (i);
        Vector_3D_impl force (atom.kit().potential_force ());

        //double fac = 0.0001 * dT * dT;
//        double factor = force.length() * limit * masses_().inverse(i);
        //double factor = force.length() * lIforce_limit_;
        double factor = force.length_2() * lIforce_limit_2;

        if (factor > 1.)
        {
            force *= 0.1 / factor;
            atom.kit().set_potential_force (force);
             
            Vector_3D_impl v (atom.kit().velocity ());
            v *= 0.9;                                       //fix
            atom.kit().set_velocity (Vector_3D_impl(v));

            if (log_it)
            {
                Text message ("!!! Force on atom ");
                message += i;
                message += " ";
                message += atom.kit(). mm_type_c_str();
                message += " scaled by ";
                message += 1. / factor;
                message += "\n";
                log () << message;

                first_bad_iteration_ = false;
            }
        }
    }
}

}//MM


