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

#include "List_interactions.h"

#include "Log.h"

#include "Create.h"
#include "Model_kit.h"
#include "anint.h"
#include "Lock.h"

#include "Geometry.h"

#include <math.h>

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


namespace MM
{
bool    List_interactions::list_on_             = true;

bool    List_interactions::solvation_           = false; 
bool    List_interactions::Sheffield_solvation_ = false; 
//double  List_interactions::Sheffield_a_         = .8; 
double  List_interactions::Sheffield_a_         = 2.; 
double  List_interactions::Sheffield_b_         = 1.;
double  List_interactions::Sheffield_RH_        = 1.2;

List_interactions::
List_interactions (Model & model)
:
    Yakub_Ronchi_interactions (model),
    pair_list_(0), is_locked_(false)
{
}


void List_interactions::
set_lock (bool on)
{
    profile ().interactions___add_force___others ();

    if (on)
    {
        cache_atoms ();
        cache_bonds ();
        cache_flag_ = true;
        //make_list ();
        //start ();
        is_locked_ = true;
    }
    else
    {
        flush_coords ();
        flush_force ();
        atom_.clear ();
        bond_.clear ();
        cache_flag_ = false;
        //finish ();
        is_locked_ = false;
    }

    profile ().lock ();
}


bool List_interactions::locked () const
{
    return is_locked_;
}

void List_interactions::
start ()
{
    Lock <List_interactions> lock(*this);

    list_valent_angles ();
    list_valent_torsions ();
    list_nonbonded_14 ();

    //if (!cache_flag_)
    //{
    //    cache_atoms ();
    //    cache_bonds ();
    //}

    list (1000000);
    //make_volumes ();//fix

    //if (!cache_flag_)
    //{
    //    atom_.clear ();
    //    bond_.clear ();
    //}
}

void List_interactions::
finish ()
{
}

void List_interactions::
make_list ()
{
    if (!list_on_)
        return;

    list (1000000);
    //return;

    if (!cache_flag_)
    {
        cache_atoms ();
        cache_bonds ();
    }

    make_volumes ();//fix

    //int n = pair_list_.size();
    //pair_list_.clear();
    //pair_list_.reserve (n);

    int i, j, count = atom_count();

    delete [] pair_list_;
    pair_list_ = new Next_atom [count*count];

    //for(j=0;  j<count;  ++j)
    //{
    //    Atom_cache & atom_j = atom_[j];

    //    miss_neighbours (atom_j);

    //    Force_field::VdWaals * vdw_ptr = 
    //        force_field_.VdW_par [atom_j.mm_type___];
    //    
    //    for (i=j+1;  i<count;  ++i)
    //    {
    //        Atom_cache & atom_i = atom_[i];
    //        
    //        if (atom_i.transparency___ == 0.)
    //            continue;

    //        Force_field::VdWaals & vdw = vdw_ptr [atom_i.mm_type___];
    //        pair_list_.push_back (Atom_pair (j, i, vdw.A, vdw.B));
    //    }
    //    
    //    process_14 (atom_j, None);
    //    restore_neighbours ();
    //}

    Force_field & ff = force_field();

    double ad_hoc_A = ff.VdW_par[Micro_object_type::any_type]
                                [Micro_object_type::any_type].A;
    double ad_hoc_B = ff.VdW_par[Micro_object_type::any_type]
                                [Micro_object_type::any_type].B;
    bool warn = !ad_hoc ();

    for(j=0;  j<count;  ++j)
    {
        Atom_cache & atom_j = atom_[j];

        miss_neighbours (atom_j);

        if (!ad_hoc() && atom_j.mm_type___ == Micro_object_type::undefined_type)
        {
            Text message = "Undefined atom type.";
            to_user().warning (message);
        }

        Force_field::VdWaals * vdw_ptr = 
            ff.VdW_par [atom_j.mm_type___];
        
        for (i=0;  i<count;  ++i)
        {
            Atom_cache & atom_i = atom_[i];
            
            double miss = 
                //atom_i.transparency___ * atom_j.transparency___ * (i != j);
                atom_i.transparency___ * atom_j.transparency___;
            if (i == j)
                miss = 0.;

            Force_field::VdWaals & vdw = vdw_ptr [atom_i.mm_type___];
            //pair_list_.push_back (Atom_pair (/*j, */i, 
            //    (float)vdw.A, (float)vdw.B, (float)miss));

            if (warn && vdw.A == ad_hoc_A && vdw.B == ad_hoc_B)
            {
                Text message = "No van der Waals parameters for types '";
                message += ::mot().itoa (atom_i.mm_type___);   
                message += "' and '";
                message += ::mot().itoa (atom_j.mm_type___);   
                message += "'.";
                to_user().warning (message);
                warn = false;
            }

            pair_list_[j*count + i] = Next_atom (/*j, */i, 
                (float)vdw.A, (float)vdw.B, (float)miss);
        }
        
        process_14 (atom_j, No);
        restore_neighbours ();
    }

    if (!cache_flag_)
    {
        atom_.clear ();
        bond_.clear ();
    }
}

void List_interactions::
list (std::deque <Atom_chunk_big> &chunk, double skin)
//list (Array <Atom_chunk> &chunk, double skin)
{
    profile ().interactions___add_force___others ();

    if (!list_on_)
        return;

    if (!cache_flag_)
        {cache_atoms ();        cache_bonds ();}

    //int count = atom_count();
    int count = atom_.size();   //fix Why?
    //int old_chunk_size = chunk.size();
    //int old_chunk_size = chunk.capacity();
    //log() << old_chunk_size << " \n";
    //chunk.clear();
    TIME ("#############");
    chunk.resize(0);
    TIME ("resize(0)    ");
    //chunk.reserve (count > old_chunk_size ? count : old_chunk_size);

    Force_field & ff      = force_field();
    double        Rm      = ff.Rm_;
    double        I_cut_2 = ff.I_cut_* ff.I_cut_;

    double  L = 0., lIL = 0., LI2 = 0.; 
    bool    box = false;
    if (model_.kit().boundary_conditions ().type_name() == "Box")
    {
        box = true;
        L = model_.kit().boundary_conditions ().box ().a_mod_rectangle();//fix
        lIL = 1. / L;
        LI2 = L / 2.;

        if (Rm > LI2)
            Rm = LI2;
    }

    double R_2 = (Rm + skin) * (Rm + skin);

    double ad_hoc_A = ff.VdW_par[Micro_object_type::any_type]
                                [Micro_object_type::any_type].A;
    double ad_hoc_B = ff.VdW_par[Micro_object_type::any_type]
                                [Micro_object_type::any_type].B;
    bool warn = !ad_hoc ();

    for(int j=0;  j<count;  ++j)
    {
        to_user().clock ();

        Atom_cache & atom_j = atom_[j];

        if (atom_j.origin___  == &Atom::none())
            break;

        miss_neighbours (atom_j);

        if (!ad_hoc() && atom_j.mm_type___ == Micro_object_type::undefined_type)
        {
            Text message = "Undefined atom type.";
            to_user().warning (message);
        }

        Force_field::VdWaals * vdw_ptr = ff.VdW_par [atom_j.mm_type___];

        double R_star = ff.R_eps[atom_j.mm_type___].R_star;
        double I_R_star_3_j = ff.I_ * R_star * R_star * R_star;
        double q2_j = atom_j.induction_charge___ * atom_j.induction_charge___;

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;

        int index = -1;

        //for (i=0;  i<count;  ++i)
        for (int k=0, i=j+1;  i<count;  ++i)
        {
            Atom_cache & atom_i = atom_[i];
            
            //index = k % big_chunk_size ();
            //if (index == 0)
            //    chunk.push_back (Atom_chunk (j, big_chunk_size()));
            
            double miss = atom_i.transparency___ * atom_j.transparency___;
            //if (i == j) miss = 0.;
            if (miss == 0.) 
            {
                //Atom_chunk & back = chunk.back();
                //back.next[index].index = count-1;
                //back.next[index].A = 0.f;
                //back.next[index].B = 0.f;
                //++k;
                continue;
            }

            double x = ax - atom_i.x___;
            double y = ay - atom_i.y___;
            double z = az - atom_i.z___;

            if (box)
            { 
                x -= L * anint (x * lIL);
                y -= L * anint (y * lIL);
                z -= L * anint (z * lIL);
            }

            double r_2 = x*x + y*y + z*z + 0.000001;
            if (r_2 > R_2)
            {
                //Atom_chunk & back = chunk.back();
                //back.next[index].index = count-1;
                //back.next[index].A = 0.f;
                //back.next[index].B = 0.f;
                //++k;
                continue;
            }

            Force_field::VdWaals & vdw = vdw_ptr [atom_i.mm_type___];

            if (warn && vdw.A == ad_hoc_A && vdw.B == ad_hoc_B)
            {
                Text message = "No van der Waals parameters for types '";
                message += ::mot().itoa (atom_i.mm_type___);   
                message += "' and '";
                message += ::mot().itoa (atom_j.mm_type___);   
                message += "'.";
                to_user().warning (message);
                warn = false;
            }

            //pair_list_[j*count + i] = Next_atom (/*j, */i, 
            //    (float)vdw.A, (float)vdw.B, (float)miss);
            ////pair_list_.push_back (Atom_pair (/*j, */i, 
            ////    (float)vdw.A, (float)vdw.B, (float)miss));
            
            //if (i % chunk_size () == 0)
            //{
                //chunk.back().next.reserve (chunk_size ());
                //chunk.back().index   = j;
            //}
            //chunk.back().next.push_back (Atom_chunk::Next());
            //chunk.back().next.back().index = i;
            //chunk.back().next.back().A     = (float)vdw.A;
            //chunk.back().next.back().B     = (float)vdw.B;
            //if (miss == 0.)
            //    chunk.back().next.back().A = -1.;

            double A = vdw.A;
            double B = vdw.B;

            double q2_i = atom_i.induction_charge___ * atom_i.induction_charge___;
            B += I_R_star_3_j * (q2_j < I_cut_2) * q2_i;

            double R_star_i = ff.R_eps[atom_i.mm_type___].R_star;
            double I_R_star_3_i = ff.I_ * R_star_i * R_star_i * R_star_i;
            B += I_R_star_3_i * (q2_i < I_cut_2) * q2_j;

            index = k % big_chunk_size ();
            if (index == 0)
                //chunk.push_back (Atom_chunk (j, big_chunk_size()));
                //chunk.push_back (Atom_chunk_big (j, big_chunk_size()));
                chunk.push_back (Atom_chunk_big (j));

            //Atom_chunk & back = chunk.back();
            Atom_chunk_big & back = chunk.back();
            back.next[index].index = i;
            back.next[index].A = (float)A;
            back.next[index].B = (float)B;
            if (miss == 0.)
                chunk.back().next[index].A = -1.;

            ++k;
        }
        for (++index;  index<big_chunk_size() && index >0;  ++index)
        {
            //Atom_chunk & back = chunk.back();
            Atom_chunk_big & back = chunk.back();
#ifndef NDEBUG
            int ind = back.next[index].index;
            float A = back.next[index].A;
            float B = back.next[index].B;
            if (!(ind == 0 && A == -1.f && B == 0.f))
                continue;
            CHECK ("void atom", ind == 0 && A == -1.f && B == 0.f);
#endif
            back.next[index].index = count-1;
            back.next[index].A = 0.f;
            back.next[index].B = 0.f;
        }
        
        process_14 (atom_j, No);
        restore_neighbours ();
    }

    if (!cache_flag_)
        {atom_.clear ();        bond_.clear ();}

    TIME ("list         ");

    profile ().interactions___add_force___make_list ();
}

void List_interactions::
list (Array      <Atom_chunk_small> &chunk, double R, 
      std::deque <Atom_chunk_big> &from)
//list (Array <Atom_chunk> &chunk, double R, 
//      Array <Atom_chunk> &from)
{
    profile ().interactions___add_force___others ();

    if (!list_on_)
        return;

    if (!cache_flag_)
        {cache_atoms ();        cache_bonds ();}

    //int old_chunk_size = chunk.size();
    //int old_chunk_size = chunk.capacity();
    //log() << old_chunk_size << " \n";
    chunk.resize(0);
    //chunk.clear();
    //chunk.reserve (atom_count() > old_chunk_size ? atom_count() : old_chunk_size);

    //Force_field & ff = force_field();

    double  L = 0., lIL = 0., LI2 = 0.; 
    bool    box = false;
    if (model_.kit().boundary_conditions ().type_name() == "Box")
    {
        box = true;
        L = model_.kit().boundary_conditions ().box ().a_mod_rectangle();//fix
        lIL = 1. / L;
        LI2 = L / 2.;
    }

    double R_2 = R * R;
    int chunk_count = from.size();

    for (int j=0;  j<chunk_count;  ++j)
    {
        to_user().clock ();

        Atom_chunk_big & chunk_j = from[j];
        int atom_j_index = chunk_j.index;
        Atom_cache & atom_j  = atom_[atom_j_index];

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;

        int count_i = chunk_j.count();

        for (int k=0,  i=0;  i<count_i;  ++i)
        {
            int atom_i_index = chunk_j.next[i].index;
            Atom_cache & atom_i = atom_[atom_i_index];
            
            double x = ax - atom_i.x___;
            double y = ay - atom_i.y___;
            double z = az - atom_i.z___;

            if (box)
            { 
                x -= L * anint (x * lIL);
                y -= L * anint (y * lIL);
                z -= L * anint (z * lIL);
            }

            double r_2 = x*x + y*y + z*z + 0.000001;
            if (r_2 > R_2)
                continue;

            double A = chunk_j.next[i].A;
            double B = chunk_j.next[i].B;

            int index = k % small_chunk_size ();
            if (index == 0)
                //chunk.push_back (Atom_chunk_small (atom_j_index, small_chunk_size()));
                chunk.push_back (Atom_chunk_small (atom_j_index));
            chunk.back().next[index].index = atom_i_index;
            chunk.back().next[index].A = (float)A;
            chunk.back().next[index].B = (float)B;

            ++k;
        }
    }

    if (!cache_flag_)
        {atom_.clear ();        bond_.clear ();}

    profile ().interactions___add_force___make_sublist ();
}

void List_interactions::
list (Array <Atom_chunk_small> &chunk, double R, 
      Array <Atom_chunk_small> &from)
//list (Array <Atom_chunk> &chunk, double R, 
//      Array <Atom_chunk> &from)
{
    profile ().interactions___add_force___others ();

    if (!list_on_)
        return;

    if (!cache_flag_)
        {cache_atoms ();        cache_bonds ();}

    //int old_chunk_size = chunk.size();
    //int old_chunk_size = chunk.capacity();
    //log() << old_chunk_size << " \n";
    chunk.resize(0);
    //chunk.clear();
    //chunk.reserve (atom_count() > old_chunk_size ? atom_count() : old_chunk_size);

    //Force_field & ff = force_field();

    double  L = 0., lIL = 0., LI2 = 0.; 
    bool    box = false;
    if (model_.kit().boundary_conditions ().type_name() == "Box")
    {
        box = true;
        L = model_.kit().boundary_conditions ().box ().a_mod_rectangle();//fix
        lIL = 1. / L;
        LI2 = L / 2.;
    }

    double R_2 = R * R;
    int chunk_count = from.size();

    for (int j=0;  j<chunk_count;  ++j)
    {
        to_user().clock ();

        Atom_chunk_small & chunk_j = from[j];
        int atom_j_index = chunk_j.index;
        Atom_cache & atom_j  = atom_[atom_j_index];

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;

        int count_i = chunk_j.count();

        for (int k=0,  i=0;  i<count_i;  ++i)
        {
            int atom_i_index = chunk_j.next[i].index;
            Atom_cache & atom_i = atom_[atom_i_index];
            
            double x = ax - atom_i.x___;
            double y = ay - atom_i.y___;
            double z = az - atom_i.z___;

            if (box)
            { 
                x -= L * anint (x * lIL);
                y -= L * anint (y * lIL);
                z -= L * anint (z * lIL);
            }

            double r_2 = x*x + y*y + z*z + 0.000001;
            if (r_2 > R_2)
                continue;

            double A = chunk_j.next[i].A;
            double B = chunk_j.next[i].B;

            int index = k % small_chunk_size ();
            if (index == 0)
                //chunk.push_back (Atom_chunk_small (atom_j_index, small_chunk_size()));
                chunk.push_back (Atom_chunk_small (atom_j_index));
            chunk.back().next[index].index = atom_i_index;
            chunk.back().next[index].A = (float)A;
            chunk.back().next[index].B = (float)B;

            ++k;
        }
    }

    if (!cache_flag_)
        {atom_.clear ();        bond_.clear ();}

    profile ().interactions___add_force___make_sublist ();
}

void List_interactions::
list_water (Atom_cache & atom_j)
{
    REQUIRE ("Exiting atom", atom_j.number___ >=0);

    water_.push_back (atom_j.number___);
}

void List_interactions::
add_water_force ()
{
    int count = water_.size();

    for (int j=0;  j<count;  ++j)
    {
        Atom_cache & atom_j = atom_[water_[j]];
        SPCf_potential (atom_j, true);
    }
}

void List_interactions::
list_valent_angles ()
{
    water_.clear ();
    angle_.clear ();

    double  K=0,  angleEq=0;

    int count = atom_count();

    for (int j=0;  j<count;  ++j)
    {
        Atom_cache & atom_j = atom_[j];

        if (atom_j.mm_type___ == water)
        {
            list_water (atom_j);
            continue;
        }

        Force_field & ff = force_field();

        int     number_of_bond = atom_j.bound_count();
        int     i, k;

        for (i=0;  i<number_of_bond;  i++)
        {
            Atom_cache const & atom_i = atom_j.bound (i);

            for (k=i+1;  k<number_of_bond;  k++)
            {
                Atom_cache const & atom_k = atom_j.bound (k);
                ff.angle (atom_i.mm_type___, atom_j.mm_type___, atom_k.mm_type___,
                    K, angleEq );

                if (ad_hoc() && ff.fail())
                {
                    K = 30.0;
                    angleEq = 125 * __Pi__ / 180;
                }

                angle_.push_back (Valent_angle_in_list (
                    atom_i.number___,
                    atom_j.number___,
                    atom_k.number___,
                    K, angleEq ));
            }
        }
    }
}

void List_interactions::
list_valent_torsions ()
{
    torsion_.clear ();

    if (&bond_group_ == &nil<Bond_group>())
        return;

    Force_field & ff = force_field();
    int         i,l, tor;
    int count = bond_count();

    for (int j=0;  j<count;  ++j)
    {
        Bond_cache & bond_j = bond_[j];
        Atom_cache & atom_j = *bond_j.first_atom___;
        Atom_cache & atom_k = *bond_j.second_atom___;

        for (i=0;  i<atom_j.bound_count();  ++i)
        {
            Atom_cache & atom_i = atom_j.bound (i);
            if (&atom_i == &atom_k) continue;

            for (l=0;  l<atom_k.bound_count();  ++l)
            {
                Atom_cache & atom_l = atom_k.bound (l);
                if (&atom_j == &atom_l) continue;
                                                               //fix may be 0             
                const Array <Force_field::Torsions> * tor_tmp; //!!!
                tor_tmp = ff.torsion (atom_i.mm_type___, 
                                      atom_j.mm_type___, 
                                      atom_k.mm_type___, 
                                      atom_l.mm_type___);
                //fix
                if (tor_tmp == 0)
                     continue;

                int tor_max = tor_tmp->size();
                for (tor=0;  tor<tor_max;  tor++)
                {
                    int    Mult  = (*tor_tmp)[tor].Mult;
                    double V_2   = (*tor_tmp)[tor].half_V;

                    if (V_2 < d_mach_eps_2)
                        continue;

                    double Phi0 = (*tor_tmp)[tor].Fi0;
                    int    n    = (*tor_tmp)[tor].N;

                    torsion_.push_back (Valent_torsion_in_list(
                        atom_i.number___,
                        atom_j.number___,
                        atom_k.number___,
                        atom_l.number___,
                        Mult, V_2, Phi0, n ));
                }
            }
        }
    }
}

void List_interactions::
list_nonbonded_14 ()
{
    nonbonded_14_.clear ();

    int count = atom_count();

    for (int j=0;  j<count;  ++j)
    {
        Atom_cache & atom_j = atom_[j];
        miss_neighbours (atom_j);
        process_14      (atom_j, List);
        restore_neighbours ();
    }
}

void List_interactions::
list_nonbonded_14 (Atom_cache & atom_1, Atom_cache & atom_2)
{
    Force_field & ff = force_field();
    double A, B;
    ff.VdW (atom_1.mm_type___, atom_2.mm_type___, A, B);

    Atom_chunk_template <double, 1>
                    chunk (atom_1.number___);
    chunk.next[0].index  = atom_2.number___;
    chunk.next[0].A      = A;
    chunk.next[0].B      = B;

    nonbonded_14_.push_back (chunk);
}

void List_interactions::
add_nonbonded_14 ()
{
    Force_field & ff = force_field();
    int count = nonbonded_14_.size();

    for (int m=0;  m<count;  ++m)
    {
        Atom_chunk_template <double, 1> const & chunk = nonbonded_14_[m];

        Atom_cache & atom_1 = atom_[chunk.index];
        Atom_cache & atom_2 = atom_[chunk.next[0].index];
        double  A = chunk.next[0].A;
        double  B = chunk.next[0].B;

        double  x = atom_1.x___ - atom_2.x___;
        double  y = atom_1.y___ - atom_2.y___;
        double  z = atom_1.z___ - atom_2.z___;

        double r_2 = x*x + y*y + z*z;

        double lIr_2 = 1. / (x*x + y*y + z*z);  // 1 / r^2
        double lIr_eps = sqrt (lIr_2) / ff.eps; // 1 / r * eps
        double lIr_6 = lIr_2 * lIr_2 * lIr_2;   // 1 / r^6

        double K  = - (ff.VdW_Scale * 
            (6 * (B - 2 * A * lIr_6) * (lIr_6 * lIr_2))
            - ff.E_Static_Scale * (atom_1.charge___ * atom_2.charge___ * EPS0 * //fix EPS0
            lIr_eps * lIr_2)
                      );

        double fx = K * x;
        double fy = K * y;
        double fz = K * z;

        atom_1.force_x___ += fx;
        atom_2.force_x___ -= fx;

        atom_1.force_y___ += fy;
        atom_2.force_y___ -= fy;

        atom_1.force_z___ += fz;
        atom_2.force_z___ -= fz;

        atom_1.F_torsion_x___ += fx;
        atom_2.F_torsion_x___ -= fx;

        atom_1.F_torsion_y___ += fy;
        atom_2.F_torsion_y___ -= fy;

        atom_1.F_torsion_z___ += fz;
        atom_2.F_torsion_z___ -= fz;

        add_virial (K*r_2, K*x*x, K*y*y, K*z*z);
    }
}

void List_interactions::
add_valent_angle_force ()
{
    int count = angle_.size();

    for (int m=0;  m<count;  ++m)
    {
        Valent_angle_in_list const & angle = angle_[m];

        Atom_cache & atom_i = atom_[angle.atom_1];
        Atom_cache & atom_j = atom_[angle.atom_2];
        Atom_cache & atom_k = atom_[angle.atom_3];
        double angleEq = angle.equilibrium_value;
        double K       = angle.force_constant;

        double  angle_minus_angleEq = 
                MM::angle (atom_i, atom_j, atom_k) -  angleEq;

        double dV_I_dTheta = 2. * K * angle_minus_angleEq;

        //-----------------------------------------------------------------
        Vector_3D_impl Ri  (atom_i.x___, atom_i.y___, atom_i.z___);
        Vector_3D_impl Rj  (atom_j.x___, atom_j.y___, atom_j.z___);
        Vector_3D_impl Rk  (atom_k.x___, atom_k.y___, atom_k.z___);
        Vector_3D_impl Rij (Rj);    Rij -= Ri; 
        Vector_3D_impl Rkj (Rj);    Rkj -= Rk;
        Vector_3D_impl M   (Rij);   M.cross (Rkj);

        //if (M.length_2() < d_mach_eps_2)
        if (M.length_2() < d_mach_eps)
            continue;

        Vector_3D_impl Fi  (Rij);   Fi.cross (M);   Fi.norm ();
        Fi *=   dV_I_dTheta / Rij.length();

        Vector_3D_impl Fk  (Rkj);   Fk.cross (M);   Fk.norm ();
        Fk *= - dV_I_dTheta / Rkj.length();

        Vector_3D_impl Fj  (Fi);    Fj += Fk;       Fj.negate ();

        atom_i.force_x___ += Fi.x();
        atom_i.force_y___ += Fi.y();
        atom_i.force_z___ += Fi.z();

        atom_j.force_x___ += Fj.x();
        atom_j.force_y___ += Fj.y();
        atom_j.force_z___ += Fj.z();

        atom_k.force_x___ += Fk.x();
        atom_k.force_y___ += Fk.y();
        atom_k.force_z___ += Fk.z();

        atom_i.F_angle_x___ += Fi.x();
        atom_i.F_angle_y___ += Fi.y();
        atom_i.F_angle_z___ += Fi.z();

        atom_j.F_angle_x___ += Fj.x();
        atom_j.F_angle_y___ += Fj.y();
        atom_j.F_angle_z___ += Fj.z();

        atom_k.F_angle_x___ += Fk.x();
        atom_k.F_angle_y___ += Fk.y();
        atom_k.F_angle_z___ += Fk.z();
    }
}

void List_interactions::
add_valent_torsion_force ()
{
    int count = torsion_.size();

    for (int m=0;  m<count;  ++m)
    {
        Valent_torsion_in_list const & torsion = torsion_[m];

        Atom_cache & atom_i = atom_[torsion.atom_1];
        Atom_cache & atom_j = atom_[torsion.atom_2];
        Atom_cache & atom_k = atom_[torsion.atom_3];
        Atom_cache & atom_l = atom_[torsion.atom_4];

        int    Mult = torsion.Mult;
        double V_2  = torsion.half_V;
        double Phi0 = torsion.Fi0;
        int    n    = torsion.N;

        double Phi = MM::torsion (atom_i, atom_j, atom_k, atom_l); 
        //Point_3D_impl  pi  (atom_i.x___, atom_i.y___, atom_i.z___);
        //Point_3D_impl  pj  (atom_j.x___, atom_j.y___, atom_j.z___);
        //Point_3D_impl  pk  (atom_k.x___, atom_k.y___, atom_k.z___);
        //Point_3D_impl  pl  (atom_l.x___, atom_l.y___, atom_l.z___);
        //double Phi =  Dihedral_angle (pi, pj, pk, pl).value(); 

        // e += V_2 * (1 + cos (n*Fi - Fi0)) / Mult;
        double dV_I_dPhi = - V_2 * n * sin (n * Phi - Phi0) / Mult;

        Vector_3D_impl  Ri  (atom_i.x___, atom_i.y___, atom_i.z___);
        Vector_3D_impl  Rj  (atom_j.x___, atom_j.y___, atom_j.z___);
        Vector_3D_impl  Rk  (atom_k.x___, atom_k.y___, atom_k.z___);
        Vector_3D_impl  Rl  (atom_l.x___, atom_l.y___, atom_l.z___);

        Vector_3D_impl  Rij (Rj);   Rij -= Ri;
        Vector_3D_impl  Rkj (Rj);   Rkj -= Rk;
        Vector_3D_impl  Rkl (Rl);   Rkl -= Rk;
        Vector_3D_impl  M   (Rij);  M.cross (Rkj);
        Vector_3D_impl  N   (Rkj);  N.cross (Rkl);
        double          kj_length_2 = Rkj.length_2 ();
        double          kj_length   = sqrt (kj_length_2);
        double          M_length_2 = M.length_2();
        double          N_length_2 = N.length_2();

        if (M_length_2 < d_mach_eps_2 || N_length_2 < d_mach_eps_2)
            continue;

        Vector_3D_impl Fi  (M);
        Fi *= - dV_I_dPhi * kj_length / M_length_2;

        Vector_3D_impl Fl  (N);
        Fl *=   dV_I_dPhi * kj_length / N_length_2;

        atom_i.force_x___ += Fi.x(); 
        atom_i.force_y___ += Fi.y(); 
        atom_i.force_z___ += Fi.z();

        atom_l.force_x___ += Fl.x(); 
        atom_l.force_y___ += Fl.y(); 
        atom_l.force_z___ += Fl.z();

        atom_i.F_torsion_x___ += Fi.x(); 
        atom_i.F_torsion_y___ += Fi.y(); 
        atom_i.F_torsion_z___ += Fi.z();

        atom_l.F_torsion_x___ += Fl.x(); 
        atom_l.F_torsion_y___ += Fl.y(); 
        atom_l.F_torsion_z___ += Fl.z();

        double Rij_dot_Rkj = Rij.dot(Rkj);  Rij_dot_Rkj /= kj_length_2;
        double Rkl_dot_Rkj = Rkl.dot(Rkj);  Rkl_dot_Rkj /= kj_length_2;

        Vector_3D_impl Fj  (Fi);    Fj.negate ();
        Vector_3D_impl Fk  (Fl);    Fk.negate ();

        Fi *= Rij_dot_Rkj;
        Fl *= Rkl_dot_Rkj;

        Fj += Fi;   Fj -= Fl;
        Fk -= Fi;   Fk += Fl;

        atom_j.force_x___ += Fj.x(); 
        atom_j.force_y___ += Fj.y(); 
        atom_j.force_z___ += Fj.z();

        atom_k.force_x___ += Fk.x(); 
        atom_k.force_y___ += Fk.y(); 
        atom_k.force_z___ += Fk.z();

        atom_j.F_torsion_x___ += Fj.x(); 
        atom_j.F_torsion_y___ += Fj.y(); 
        atom_j.F_torsion_z___ += Fj.z();

        atom_k.F_torsion_x___ += Fk.x(); 
        atom_k.F_torsion_y___ += Fk.y(); 
        atom_k.F_torsion_z___ += Fk.z();
    }
}

void List_interactions::
cache_coords ()
{
    REQUIRE ("Atoms cached.", locked());

    int i, count = atom_group_.atom_count();
    
    for (i=0;  i<count;  ++i)
    {
        Atom & current = atom_group_.atom (i);
        Atom_cache & cache = atom_[i];
        cache.x___ = current.x();
        cache.y___ = current.y();
        cache.z___ = current.z();
    }
}


double List_interactions::
potential ()
{
    //Lock

    int i, count;
    last_potential_  = 0.;

    Van_der_Waals_potential_    = 0.;
    electrostatic_potential_    = 0.;
    bond_potential_             = 0.;
    angle_potential_            = 0.;
    torsion_potential_          = 0.;
    improper_torsion_potential_ = 0.;
    hydrophobic_potential_      = 0.;
    surface_area_potential_     = 0.;
    generalized_Born_potential_ = 0.;
    induction_potential_        = 0.;
    water_potential_            = 0.;
    intermolecular_Van_der_Waals_potential_ = 0.;
    intermolecular_electrostatic_potential_ = 0.;
    intermolecular_induction_potential_     = 0.;

    try
    {
        if (!cache_flag_)
        {
            cache_atoms ();
            cache_bonds ();

            if (!locked())
            {
                cache_flag_ = true; //fix !!! always
                list (1000000);
                //make_list   ();
                cache_flag_ = false;
            }
        }
    }
    catch(...)
    {
        to_user().fatal_error ("List_interactions:\n"
                               "Can not make atoms pair list.\n"
                               "Possible reason is a lack of memory.");
        return 0;
    }

    try
    {
        full_nonbonded_potential ();

        count = atom_count();
        for (i=0;  i<count;  ++i)
        {
            Atom_cache & current_atom = atom_[i];
            atom_potential (current_atom);
        }

        if (&bond_group_ != &nil<Bond_group>())
        {
            count = bond_count();
            for (i=0;  i<count;  ++i)
            {
                Bond_cache & current_bond = bond_[i];             
                bond_potential (current_bond);
            }
        }
    }
    catch(...)
    {
        to_user().fatal_error ("List_interactions: Unknown exception.");
    }

    if (!cache_flag_)
    {
        //flush_born  ();//fix down
        atom_.clear ();
        bond_.clear ();
    }

    return last_potential_;
}

void List_interactions::
add_force ()
{
    Timer timer;
    Text  time;

    int i, count;

    //profile().interactions___add_force___start ();

    if (!cache_flag_)
    {
        cache_atoms ();
        cache_bonds ();

        if (!locked())
        {
            cache_flag_ = true; //fix !!! always
            //make_list   ();
            list (1000000);
            cache_flag_ = false;

            //Conventional_interactions::add_force ();
            //return;
        }
    }

    //profile().interactions___add_force___make_list ();

    time += "   cache ";
    time += timer.elapsed();    timer.restart();

    try
    {
        zero_forces ();      //fix remove?
        set_virial (0,0,0,0);

        calc_nonbonded_force ();

        //profile().interactions___add_force___calc_nonbonded_force ();


        time += "   nb ";
        time += timer.elapsed();    timer.restart();

        profile().interactions___add_force___others ();
        add_water_force ();
        profile().interactions___add_force___add_water_force ();

        profile().interactions___add_force___others ();
        add_valent_angle_force ();

        count = atom_count ();
        for (i=0;  i<count;  ++i)
            add_atom_force (atom_[i]);
        profile().interactions___add_force___add_atom_force ();

        time += "   atom ";
        time += timer.elapsed();    timer.restart();

        add_valent_torsion_force ();
        add_nonbonded_14 ();

        if (&bond_group_ != &nil<Bond_group>())
        {
            count = bond_count ();

            for (i=0;  i<count;  ++i)
                add_bond_force (bond_[i]);
        }
        profile().interactions___add_force___add_bond_force ();

        time += "   bond ";
        time += timer.elapsed();    timer.restart();

        add_full_nonbonded_force ();

        //profile().interactions___add_force___add_full_nonbonded_force ();

        time += "   add nb ";
        time += timer.elapsed();    timer.restart();
    }
    catch(...)
    {
        to_user().fatal_error ("List_interactions: Unknown exception.");
    }

    if (!cache_flag_)
    {
        flush_force ();
        //flush_born  ();//fix down
        atom_.clear ();
        bond_.clear ();
    }

    //profile().List_interactions___add_force___flush_force ();

    time += "   flush ";
    time += timer.elapsed();    timer.restart();
    //to_user().status (time);
    to_user().clock ();
}

void List_interactions::
add_short_force ()
{
    int i, count;
    TIME ("add_short_force              ### ");

    try
    {
        //profile().interactions___add_force___others ();
        zero_forces ();      //fix remove?

        profile().interactions___add_force___others ();
        add_water_force ();
        profile().interactions___add_force___add_water_force ();

        count = atom_count ();

        //profile().interactions___add_force___others ();
    TIME ("zero_forces                      ");
        add_valent_angle_force ();

        for (i=0;  i<count;  ++i)
            add_atom_force (atom_[i]);
    TIME ("add_atom_force                   ");
        profile().interactions___add_force___add_atom_force ();

        add_valent_torsion_force ();
        add_nonbonded_14 ();

        if (&bond_group_ != &nil<Bond_group>())
        {
            count = bond_count ();

            for (i=0;  i<count;  ++i)
                add_bond_force (bond_[i]);
        }
        profile().interactions___add_force___add_bond_force ();
    }
    catch(...)
    {
        to_user().fatal_error ("List_interactions: Unknown exception.");
    }
    TIME ("add_short_force              ___ ");
}

void List_interactions::
full_nonbonded_potential ()
{
    Timer timer;

    if (Sheffield_solvation_)
        Sheffield_potential ();

    if (!is_on())
        return;

    Force_field & ff = force_field();
    double Rm      = ff.Rm_;
    double I_cut_2 = ff.I_cut_* ff.I_cut_;

    double  L = 0., lIL = 0., LI2 = 0.;
    bool    box = false;
    if (model_.kit().boundary_conditions ().type_name() == "Box")
    {
        box = true;
        L = model_.kit().boundary_conditions ().box ().a_mod_rectangle();//fix
        lIL = 1. / L;
        LI2 = L / 2.;

        if (Rm > LI2)
            Rm = LI2;
    }

    const double Rm_2    = Rm * Rm;
    const double lIRm    = 1. / Rm;

    double  lIepsilon = Constant::lI4PiEps0 / epsilon();
    //double  A = 0., B = 0.;
    double  x, y, z;
    double  r_2, lIr, lIr_6, lIr_eps;
    double  P_VdW = 0., P_el = 0., Yakub_Ronchi = 0.;
    
    //fix I
    double  Ind = 0.;

    //int count = pair_list_.size();

    //for (int i=0;  i<count;  ++i)
    //{
    //    Atom_pair const & pair = pair_list_[i];
    //    Atom_cache & atom1 = atom_[pair.atom1_N___];
    //    Atom_cache & atom2 = atom_[pair.atom2_N___];

    //    A = pair.A___;
    //    B = pair.B___;

    //    x = atom1.x___ - atom2.x___;
    //    y = atom1.y___ - atom2.y___;
    //    z = atom1.z___ - atom2.z___;

    //    r_2     = x*x + y*y + z*z;
    //    lIr     = 1.0 / sqrt (r_2);
    //    lIr_eps = lIr * lIepsilon;
    //    lIr_6   = lIr * lIr;    lIr_6 *= lIr_6 * lIr_6;

    //    P_VdW += (A * lIr_6 - B ) * lIr_6;
    //    P_el  += atom1.charge___ * atom2.charge___ * lIr_eps; //fix
    //}

    //int count = atom_count();
    int chunk_count = chunk_.size();

    //for (int j=0;  j<count;  ++j)
    for (int j=0;  j<chunk_count;  ++j)
    {
        to_user().clock ();

        //Atom_cache & atom_j = atom_[j];
        Atom_chunk_big & chunk_j = chunk_[j];
        int atom_j_index = chunk_j.index;
        Atom_cache & atom_j  = atom_[atom_j_index];
        
        //fix I
        double R_star = ff.R_eps[atom_j.mm_type___].R_star;
        double I_R_star_3_j = ff.I_ * R_star * R_star * R_star;
        double q2_j = atom_j.induction_charge___ * atom_j.induction_charge___;

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;

        //int count_i = j*count + count;
        int count_i = chunk_j.count();

        for (int i=0;  i<count_i;  ++i)
        //for (int i=j*count + j +1 ;  i<count_i;  ++i)
        {
            //Next_atom const & pair = pair_list_[i];
            //Atom_cache & atom_i = atom_[pair.atom2_N___];
            //A = pair.A___;
            //B = pair.B___;
            Atom_chunk_big::Next & next = chunk_j.next[i];
            int atom_i_index = next.index;
            double A = next.A;
            double B = next.B;
            if (A < 0) continue;
            Atom_cache & atom_i = atom_[atom_i_index];

            x = ax - atom_i.x___;
            y = ay - atom_i.y___;
            z = az - atom_i.z___;

            if (box)
            {
                x -= L * anint (x * lIL);
                y -= L * anint (y * lIL);
                z -= L * anint (z * lIL);
            }

            r_2     = x*x + y*y + z*z + 0.000001;
            lIr     = /*pair.miss___ * */1.0 / sqrt (r_2);
            lIr_eps = lIr * lIepsilon;
            lIr_6   = lIr * lIr;    lIr_6 *= lIr_6 * lIr_6;

            double rIRm = r_2 * lIr * lIRm;
            Yakub_Ronchi   = 1. + .5 * rIRm * (rIRm * rIRm - 3.); 

            double current_P_VdW = (A * lIr_6 - B ) * lIr_6;
            P_VdW += current_P_VdW;

            double current_P_el = (r_2 < Rm_2) * Yakub_Ronchi *
                atom_j.charge___ * atom_i.charge___ * lIr_eps; //fix
            P_el += current_P_el;

            //fix I
            //double current_Ind = -I_R_star_3 * atom2.charge___ * atom2.charge___ * lIr_6;
            double q2_i = atom_i.induction_charge___ * atom_i.induction_charge___;
            double current_Ind = 
                -I_R_star_3_j * (q2_j < I_cut_2) * q2_i * lIr_6;

            double R_star_i = ff.R_eps[atom_i.mm_type___].R_star;
            double I_R_star_3_i = ff.I_ * R_star_i * R_star_i * R_star_i;
            //current_Ind += -I_R_star_3_i * atom1.charge___ * atom1.charge___ * lIr_6;
            current_Ind += 
                -I_R_star_3_i * (q2_i < I_cut_2) * q2_j * lIr_6;

            Ind += current_Ind;

            if (atom_j.mol___ != atom_i.mol___)
            {
                intermolecular_Van_der_Waals_potential_ += current_P_VdW;
                intermolecular_electrostatic_potential_ += current_P_el;
                intermolecular_induction_potential_     += current_Ind;
            }
        }
    }

    Van_der_Waals_potential_ += P_VdW;
    electrostatic_potential_ += P_el;
    //hydrophobic_potential_   += P_hph;
    induction_potential_     += Ind;

    //fix I
    last_potential_ += P_el + P_VdW /*+ P_hph*/ + Ind;

    //log() << "\n Nonbonded " << timer.elapsed();
}

void List_interactions::
add_full_nonbonded_force ()
{
    profile().interactions___add_force___others ();

    if (Sheffield_solvation_)
        add_full_Sheffield_force ();

    if (!is_on())
        return;

    Force_field & ff = force_field();
    double Rm      = ff.Rm_;
    //double I_cut_2 = ff.I_cut_* ff.I_cut_; //fix here and in add_long_force (Array <Atom_chunk> & chunk)

    double  L = 0., lIL = 0., LI2 = 0.; 
    bool    box = false;
    if (model_.kit().boundary_conditions ().type_name() == "Box")
    {
        box = true;
        L = model_.kit().boundary_conditions ().box ().a_mod_rectangle();//fix
        lIL = 1. / L;
        LI2 = L / 2.;

        if (Rm > LI2)
            Rm = LI2;
    }

    const double Rm_2    = Rm * Rm;
    const double lIRm_3  = 1. / (Rm*Rm*Rm);

    //double  lIepsilon = EPS0 / epsilon();
    double  lIepsilon = Constant::lI4PiEps0 / epsilon();

    //for (int j=0;  j<count;  ++j)
    //{

    //            //if      (x >  LI2) x -= L;
    //            //else if (x < -LI2) x += L;
    //            //if      (y >  LI2) y -= L;
    //            //else if (y < -LI2) y += L;
    //            //if      (z >  LI2) z -= L;
    //            //else if (z < -LI2) z += L;

    //}

    int chunk_count = chunk_.size();

    for (int j=0;  j<chunk_count;  ++j)
    {
        to_user().clock ();

        Atom_chunk_big & chunk_j = chunk_[j];
        int atom_j_index = chunk_j.index;
        Atom_cache & atom_j  = atom_[atom_j_index];

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;
        //double Qj = atom_j.charge___;
        double QjIepsilon = atom_j.charge___* lIepsilon;

        int count_i = chunk_j.count();

        for (int i=0;  i<count_i;  ++i)
        {
            Atom_chunk_big::Next & next = chunk_j.next[i];
            int atom_i_index = next.index;
            double A = next.A;
            double B = next.B;
            if (A < 0) continue;
            Atom_cache & atom_i = atom_[atom_i_index];

            double x = ax - atom_i.x___;
            double y = ay - atom_i.y___;
            double z = az - atom_i.z___;

            if (box)
            { 
                x -= L * anint (x * lIL);
                y -= L * anint (y * lIL);
                z -= L * anint (z * lIL);
            }


            double r_2     = x*x + y*y + z*z + 0.000001;
            double lIr     = (1.0 / sqrt (r_2));
            double lIr_2   = lIr * lIr;    
            double lIr_6   = lIr_2 * lIr_2 * lIr_2;

            double Kf_Yakub_Ronchi = lIRm_3 - lIr_2 * lIr; 

            //double q = - Qj /*atom_j.charge___*/ * atom_i.charge___ * lIepsilon
            double q = - QjIepsilon * atom_i.charge___ 
                * Kf_Yakub_Ronchi * (r_2 < Rm_2);

            double VdW = - (6 * (B - 2. * A * lIr_6) * lIr_6) * lIr_2;

            double K = q + VdW /*+ Ind*/;

            double fx = K * x;
            double fy = K * y;
            double fz = K * z;

            atom_j.force_x___ += fx;
            atom_i.force_x___ -= fx;

            atom_j.force_y___ += fy;
            atom_i.force_y___ -= fy;
    
            atom_j.force_z___ += fz;
            atom_i.force_z___ -= fz;

            atom_j.F_nonbonded_x___ += fx;
            atom_i.F_nonbonded_x___ -= fx;

            atom_j.F_nonbonded_y___ += fy;
            atom_i.F_nonbonded_y___ -= fy;

            atom_j.F_nonbonded_z___ += fz;
            atom_i.F_nonbonded_z___ -= fz;

            add_virial (K*r_2, K*x*x, K*y*y, K*z*z);
        }
    }

    profile().interactions___add_force___add_full_nonbonded_force ();
}

void List_interactions::
add_long_force (std::deque <Atom_chunk_big> & chunk)
//add_long_force (Array <Atom_chunk> & chunk)
{
    profile().interactions___add_force___others ();

    zero_forces ();      //fix remove?
    //set_virial (0,0,0,0);

    if (Sheffield_solvation_)
        add_full_Sheffield_force ();

    if (!is_on())
        return;

    Force_field & ff      = force_field();
    double        Rm      = ff.Rm_;
    //double        I_cut_2 = ff.I_cut_* ff.I_cut_;

    double  L = 0., lIL = 0., LI2 = 0.; 
    bool    box = false;
    if (model_.kit().boundary_conditions ().type_name() == "Box")
    {
        box = true;
        L = model_.kit().boundary_conditions ().box ().a_mod_rectangle();//fix
        lIL = 1. / L;
        LI2 = L / 2.;

        if (Rm > LI2)
            Rm = LI2;
    }

    const double Rm_2   = Rm * Rm;
    const double lIRm_3 = 1. / (Rm*Rm*Rm);

    double  lIepsilon = Constant::lI4PiEps0 / epsilon();

    //__m128d mm_1    = _mm_set_pd  (1., 1.);
    //__m128d mm_0    = _mm_set_pd  (0.000001, 0.000001);

    int chunk_count = chunk.size();

    for (int j=0;  j<chunk_count;  ++j)
    {
        to_user().clock ();

        Atom_chunk_big & chunk_j = chunk[j];
        int atom_j_index = chunk_j.index;
        Atom_cache & atom_j  = atom_[atom_j_index];

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;
        //double Qj = atom_j.charge___;
        double QjIepsilon = atom_j.charge___* lIepsilon;

        int count_i = chunk_j.count();

        //for (int i=0;  i<count_i;  ++i)
        for (int i=0;  i<count_i;  i+=2)
        {
            //0.645 

            //_mm_prefetch ((char *)
                //&atom_[0] + chunk_j.next[i+2].index
                //, _MM_HINT_NTA); //0.680
                //, _MM_HINT_T0); //0.683
                //, _MM_HINT_T1); //0.682
                //, _MM_HINT_T2); //0.682

                //&chunk_j.next[0] + 2
                //, _MM_HINT_NTA); //
                //, _MM_HINT_T0); //
                //, _MM_HINT_T1); //0.651
                //, _MM_HINT_T2); //0.649


            Atom_chunk_big::Next & next = chunk_j.next[i];
            int atom_i_index = next.index;
            double A = next.A;
            double B = next.B;
            bool yes = A >= 0;
            //if (A < 0) 
            //    continue;
            CHECK ("No misses atoms", A>=0);
            Atom_cache & atom_i = atom_[atom_i_index];

            Atom_chunk_big::Next & next2 = chunk_j.next[i+1];
            int atom_i_index2 = next2.index;
            double A2 = next2.A;
            double B2 = next2.B;
            bool yes2 = A2 >= 0;
            //if (A2 < 0) 
            //    continue;
            CHECK ("No misses atoms", A2>=0);
            Atom_cache & atom_i2 = atom_[atom_i_index2];

            double x = ax - atom_i.x___;
            double y = ay - atom_i.y___;
            double z = az - atom_i.z___;

            double x2 = ax - atom_i2.x___;
            double y2 = ay - atom_i2.y___;
            double z2 = az - atom_i2.z___;

            if (box)
            { 
                x -= L * anint (x * lIL);
                y -= L * anint (y * lIL);
                z -= L * anint (z * lIL);

                x2 -= L * anint (x2 * lIL);
                y2 -= L * anint (y2 * lIL);
                z2 -= L * anint (z2 * lIL);
            }

            double r_2      = x*x + y*y + z*z + 0.000001;
            double r_22     = x2*x2 + y2*y2 + z2*z2 + 0.000001;

            double lIr      = (1.0 / sqrt (r_2));
            double lIr2     = (1.0 / sqrt (r_22));

/*
            //0.607 0.592
            __m128d mm_x    = _mm_set_pd  (x2, x);
            __m128d mm_y    = _mm_set_pd  (y2, y);
            __m128d mm_z    = _mm_set_pd  (z2, z);
            //__m128d mm_1    = _mm_set_pd  (1., 1.);

            __m128d mm_x2   = _mm_mul_pd (mm_x, mm_x);
            __m128d mm_y2   = _mm_mul_pd (mm_y, mm_y);
            __m128d mm_z2   = _mm_mul_pd (mm_z, mm_z);

            __m128d mm_r2   = _mm_add_pd (mm_x2, mm_y2);
                    mm_r2   = _mm_add_pd (mm_r2, mm_z2);
                    mm_r2   = _mm_add_pd (mm_r2, mm_0 );

            //__m128d mm_r_2  = _mm_set_pd  (r_22, r_2);
            __m128d mm_sqrt = _mm_sqrt_pd (mm_r2);
            __m128d mm_lIr  = _mm_div_pd  (mm_1, mm_sqrt);
            __declspec(align(16)) double r2[2];
            __declspec(align(16)) double p [2];
            _mm_store_pd (r2, mm_r2);
            _mm_store_pd (p,  mm_lIr);
            double r_2  = r2[0];
            double r_22 = r2[1];
            double lIr  = p[0];
            double lIr2 = p[1];
//*/
            double lIr_2    = lIr * lIr;    
            double lIr_22   = lIr2 * lIr2; 

            double lIr_6    = lIr_2 * lIr_2 * lIr_2;
            double lIr_62   = lIr_22 * lIr_22 * lIr_22;

            double Kf_Yakub_Ronchi = lIRm_3 - lIr_2 * lIr; 

            //double q = - Qj /*atom_j.charge___*/ * atom_i.charge___ * lIepsilon
            double q = - QjIepsilon * atom_i.charge___ 
                * Kf_Yakub_Ronchi * (r_2 < Rm_2);

            double VdW = - (6 * (B - 2. * A * lIr_6) * lIr_6) * lIr_2;

            double K = q + VdW /*+ Ind*/;



            double Kf_Yakub_Ronchi2 = lIRm_3 - lIr_22 * lIr2; 

            //double q = - Qj /*atom_j.charge___*/ * atom_i.charge___ * lIepsilon
            double q2 = - QjIepsilon * atom_i2.charge___ 
                * Kf_Yakub_Ronchi2 * (r_22 < Rm_2);

            double VdW2 = - (6 * (B2 - 2. * A2 * lIr_62) * lIr_62) * lIr_22;

            double K2 = q2 + VdW2 /*+ Ind*/;

            K  *= yes;
            K2 *= yes2;

            double fx = K * x;
            double fy = K * y;
            double fz = K * z;

            double fx2 = K2 * x2;
            double fy2 = K2 * y2;
            double fz2 = K2 * z2;

            atom_j.F_nonbonded_x___ += fx;
            atom_i.F_nonbonded_x___ -= fx;

            atom_j.F_nonbonded_y___ += fy;
            atom_i.F_nonbonded_y___ -= fy;

            atom_j.F_nonbonded_z___ += fz;
            atom_i.F_nonbonded_z___ -= fz;


            atom_j. F_nonbonded_x___ += fx2;
            atom_i2.F_nonbonded_x___ -= fx2;

            atom_j. F_nonbonded_y___ += fy2;
            atom_i2.F_nonbonded_y___ -= fy2;

            atom_j. F_nonbonded_z___ += fz2;
            atom_i2.F_nonbonded_z___ -= fz2;

            add_virial (K*r_2, K*x*x, K*y*y, K*z*z);

            add_virial (K2*r_22, K2*x2*x2, K2*y2*y2, K2*z2*z2);
        }
    }

    profile().interactions___add_force___add_long_force ();
}

void List_interactions::
add_long_force (double R, double Switching, Array <Atom_chunk_small> & chunk)
//add_long_force (double R, double Switching, Array <Atom_chunk_small> & chunk)
{
    profile().interactions___add_force___others ();

    zero_forces ();      
    //set_virial (0,0,0,0);//fix remove?

    if (!is_on())
        return;

    Force_field & ff = force_field();
    double Rm      = ff.Rm_;

    double  L = 0., lIL = 0., LI2 = 0.; 
    bool    box = false;
    if (model_.kit().boundary_conditions ().type_name() == "Box")
    {
        box = true;
        L = model_.kit().boundary_conditions ().box ().a_mod_rectangle();//fix
        lIL = 1. / L;
        LI2 = L / 2.;

        if (Rm > LI2)
            Rm = LI2;

        //if (Rm > R)
        //    Rm = R;
    }


    double       R_2 = R * R;
    const double Rm_2    = Rm * Rm;
    const double lIRm_3  = 1. / (Rm*Rm*Rm);

    // for Switching function
    double s_off = R;
    double s_off_2 = R_2;
    double s_on  = s_off - Switching;
    double s_on_2  = s_on*s_on;

    double  lIepsilon = Constant::lI4PiEps0 / epsilon();

    int chunk_count = chunk.size();

    for (int j=0;  j<chunk_count;  ++j)
    {
        Atom_chunk_small & chunk_j = chunk[j];
        int atom_j_index = chunk_j.index;
        Atom_cache & atom_j  = atom_[atom_j_index];

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;

        int count_i = chunk_j.count();

        //for (int i=0;  i<count_i;  ++i)
        for (int i=0;  i<count_i;  i+=2)
        {
            Atom_chunk_small::Next & next  = chunk_j.next[i];
            Atom_chunk_small::Next & next2 = chunk_j.next[i+1];

            int atom_i_index  = next.index;
            int atom_i_index2 = next2.index;

            double A = next.A;
            double B = next.B;
            bool yes = A >= 0;
            //if (A < 0) continue;
            Atom_cache & atom_i = atom_[atom_i_index];

            double A2 = next2.A;
            double B2 = next2.B;
            bool yes2 = A2 >= 0;
            //if (A2 < 0) continue;
            Atom_cache & atom_i2 = atom_[atom_i_index2];

            double x = ax - atom_i.x___;
            double y = ay - atom_i.y___;
            double z = az - atom_i.z___;

            double x2 = ax - atom_i2.x___;
            double y2 = ay - atom_i2.y___;
            double z2 = az - atom_i2.z___;

            if (box)
            { 
                x -= L * anint (x * lIL);
                y -= L * anint (y * lIL);
                z -= L * anint (z * lIL);

                x2 -= L * anint (x2 * lIL);
                y2 -= L * anint (y2 * lIL);
                z2 -= L * anint (z2 * lIL);
            }

            double r_2     = x*x + y*y + z*z + 0.000001;
            //if (r_2 > R_2)
            //    continue;
            double lIr     = (1.0 / sqrt (r_2));
            double lIr_2   = lIr * lIr;    
            double lIr_6   = lIr_2 * lIr_2 * lIr_2;

            double r_22     = x2*x2 + y2*y2 + z2*z2 + 0.000001;
            //if (r_2 > R_2)
            //    continue;
            double lIr2     = (1.0 / sqrt (r_22));
            double lIr_22   = lIr2 * lIr2;    
            double lIr_62   = lIr_22 * lIr_22 * lIr_22;

            double Kf_Yakub_Ronchi = lIRm_3 - lIr_2 * lIr; 

            double q = - atom_j.charge___ * atom_i.charge___ * lIepsilon
                * Kf_Yakub_Ronchi * (r_2 < Rm_2);

            double VdW = - (6 * (B - 2. * A * lIr_6) * lIr_6) * lIr_2;

            double K = q + VdW /*+ Ind*/;



            double Kf_Yakub_Ronchi2 = lIRm_3 - lIr_22 * lIr2; 

            double q2 = - atom_j.charge___ * atom_i2.charge___ * lIepsilon
                * Kf_Yakub_Ronchi2 * (r_22 < Rm_2);

            double VdW2 = - (6 * (B2 - 2. * A2 * lIr_62) * lIr_62) * lIr_22;

            double K2 = q2 + VdW2 /*+ Ind*/;

            // Switching function

            double eta = (r_2 - s_on_2) / (s_off_2 - s_on_2);
            double S;
            if (r_2 > s_off_2)
                S = 0.;
            else if (r_2 <= s_on_2)
                S = 1.;
            else
                S = 1. + eta*eta * (2*eta - 3) ;
            K *= S * yes;

            // Switching function
            double eta2 = (r_22 - s_on_2) / (s_off_2 - s_on_2);
            double S2;
            if (r_22 > s_off_2)
                S2 = 0.;
            else if (r_22 <= s_on_2)
                S2 = 1.;
            else
                S2 = 1. + eta2*eta2 * (2*eta2 - 3) ;
            K2 *= S2 * yes2;

            double fx = K * x;
            double fy = K * y;
            double fz = K * z;

            double fx2 = K2 * x2;
            double fy2 = K2 * y2;
            double fz2 = K2 * z2;

            atom_j.F_nonbonded_x___ += fx;
            atom_i.F_nonbonded_x___ -= fx;

            atom_j.F_nonbonded_y___ += fy;
            atom_i.F_nonbonded_y___ -= fy;

            atom_j.F_nonbonded_z___ += fz;
            atom_i.F_nonbonded_z___ -= fz;


            atom_j.F_nonbonded_x___  += fx2;
            atom_i2.F_nonbonded_x___ -= fx2;

            atom_j.F_nonbonded_y___  += fy2;
            atom_i2.F_nonbonded_y___ -= fy2;

            atom_j.F_nonbonded_z___  += fz2;
            atom_i2.F_nonbonded_z___ -= fz2;

            //add_virial (K*r_2, K*x*x, K*y*y, K*z*z);
        }
    }
    profile().interactions___add_force___add_switching_force ();
}



void List_interactions::
atom_potential (Atom_cache & atom)
{
    if (!is_on())
    {
        Conventional_interactions::atom_potential (atom);
        return;
    }

    miss_neighbours (atom);
    process_14      (atom, Potential);//fix ???
    restore_neighbours ();


    angle_potential            (atom);
    improper_torsion_potential (atom);
}

void List_interactions::
add_atom_force (Atom_cache & atom)
{
    //profile().interactions___add_force___others ();

    if (!is_on())
    {
        Conventional_interactions::add_atom_force (atom);
        return;
    }

    //Force_field & ff = force_field();

    ////TIME ("add_atom_force              ### ");
    //miss_neighbours (atom);
    ////TIME ("miss_neighbours                 ");
    //process_14      (atom, Force);//fix ???
    ////TIME ("process_14                      ");
    //restore_neighbours ();
    ////TIME ("restore_neighbours              ");

    ////add_angle_force            (atom, ff);
    //TIME ("add_angle_force                 ");
    add_improper_torsion_force (atom);
    //TIME ("add_improper_torsion_force      ");

    //profile().interactions___add_force___add_atom_force ();
}

void List_interactions::
add_bond_force (Bond_cache & bond)
{
    add_bond_stretch_force (bond);
    ////add_torsion_force      (bond);
}

void List_interactions::
Sheffield_potential ()
{
    int     j;
    double          e           = 0.;
    double          e_sum       = 0.;
    int             count  = atom_count();
    //int             chunk_count = chunk_.size();
    double          k = -(Constant::lI4PiEps0 / 2.) 
        * (1. - 1. / Constant::Water_dielectric_permittivity);
    Force_field &   ff = force_field();
    
    //Atom_cache * las    t_j = 0;

    //// self
    //for (j=0;  j<count;  ++j)
    //{
    //    Atom_cache & atom_j = atom_[j];
    //    double k_charge_j = k * atom_j.charge___;
    //    double R_star_j   = ff.R_eps[atom_j.mm_type___].R_star;
    //    double r_Born_2 = Sheffield_a_ * R_star_j * R_star_j + 0.000001;
    //    e = (k_charge_j * atom_j.charge___) / sqrt (r_Born_2);
    //    e_sum += e;
    //}

    //for (j=0;  j<chunk_count;  ++j)
    for(j=0;  j<count;  ++j)
    {
        to_user().clock ();

        //Atom_chunk & chunk_j = chunk_[j];
        //int atom_j_index     = chunk_j.index;
        //Atom_cache & atom_j  = atom_[atom_j_index];
        Atom_cache & atom_j = atom_[j];

        double k_charge_j = k * atom_j.charge___;
        double R_star_j   = ff.R_eps[atom_j.mm_type___].R_star;
        R_star_j = R_star_j < Sheffield_RH_ ? Sheffield_RH_ : R_star_j;

        // self
        //if (last_j != & atom_j)
        {
            double r_Born_2 = Sheffield_a_ * R_star_j * R_star_j + 0.000001;
            e = (k_charge_j * atom_j.charge___) / sqrt (r_Born_2);
            e_sum += e;
        }
        //last_j = & atom_j;

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;

        //int count_i = chunk_j.count();

        //for (int i=0;  i<count_i;  ++i)
        for (int i=j+1;  i<count;  ++i)
        {
            //Atom_chunk::Next & next = chunk_j.next[i];
            //int atom_i_index = next.index;
            ////double A = next.A;
            ////if (A < 0) continue;
            //Atom_cache & atom_i = atom_[atom_i_index];
            Atom_cache & atom_i = atom_[i];

            double R_star_i = ff.R_eps[atom_i.mm_type___].R_star;
            R_star_i = R_star_i < Sheffield_RH_ ? Sheffield_RH_ : R_star_i;
    
            double x = ax - atom_i.x___;
            double y = ay - atom_i.y___;
            double z = az - atom_i.z___;

            double r_2 = x*x + y*y + z*z + 0.000001;
            double r_Born_2  = Sheffield_a_ * R_star_j * R_star_i + 0.000001;
            double divisor_2 = r_Born_2 + /*Sheffield_b_ **/ r_2;
            double divisor   = sqrt (divisor_2);
            e = (k_charge_j * atom_i.charge___) / divisor;
            e_sum += 2 * e;
        }
    }

    generalized_Born_potential_ += e_sum;
    last_potential_ += e_sum;
}

void List_interactions::
add_full_Sheffield_force ()
{
    double          e           = 0.;
    int             count  = atom_count();
    //int             chunk_count = chunk_.size();
    double          k = -(Constant::lI4PiEps0 / 2.) 
        * (1. - 1. / Constant::Water_dielectric_permittivity);
    Force_field &   ff = force_field();
    
    //Atom_cache * last_j = 0;

    //for (int j=0;  j<chunk_count;  ++j)
    for(int j=0;  j<count;  ++j)
    {
        to_user().clock ();

        //Atom_chunk & chunk_j = chunk_[j];
        //int atom_j_index     = chunk_j.index;
        //Atom_cache & atom_j  = atom_[atom_j_index];
        Atom_cache & atom_j = atom_[j];

        double k_charge_j = k * atom_j.charge___;
        double R_star_j   = ff.R_eps[atom_j.mm_type___].R_star;
        R_star_j = R_star_j < Sheffield_RH_ ? Sheffield_RH_ : R_star_j;

        double  ax = atom_j.x___;
        double  ay = atom_j.y___;
        double  az = atom_j.z___;

        //int count_i = chunk_j.count();

        //for (int i=0;  i<count_i;  ++i)
        for (int i=j+1;  i<count;  ++i)
        {
            //Atom_chunk::Next & next = chunk_j.next[i];
            //int atom_i_index = next.index;
            ////double A = next.A;
            ////if (A < 0) continue;
            //Atom_cache & atom_i = atom_[atom_i_index];
            Atom_cache & atom_i = atom_[i];

            double R_star_i = ff.R_eps[atom_i.mm_type___].R_star;
            R_star_i = R_star_i < Sheffield_RH_ ? Sheffield_RH_ : R_star_i;
    
            double x = ax - atom_i.x___;
            double y = ay - atom_i.y___;
            double z = az - atom_i.z___;

            double r_2 = x*x + y*y + z*z + 0.000001;
            //double r   = sqrt (r_2);
            double r_Born_2  = Sheffield_a_ * R_star_j * R_star_i + 0.000001;
            double lIdivisor_2  = 1. / (r_Born_2 + /*Sheffield_b_ **/ r_2);
            double lIdivisor  = sqrt (lIdivisor_2);
            e = (k_charge_j * atom_i.charge___) * lIdivisor;
            e *= 2.; 
            double de = -e * /*r **/ lIdivisor_2;

            //de /= r;
            double K = -de;//fix

            double fx = K * x;
            double fy = K * y;
            double fz = K * z;

            atom_j.force_x___ += fx;
            atom_i.force_x___ -= fx;

            atom_j.force_y___ += fy;
            atom_i.force_y___ -= fy;

            atom_j.force_z___ += fz;
            atom_i.force_z___ -= fz;


            atom_j.F_nonbonded_x___ += fx;
            atom_i.F_nonbonded_x___ -= fx;
                
            atom_j.F_nonbonded_y___ += fy;
            atom_i.F_nonbonded_y___ -= fy;
                
            atom_j.F_nonbonded_z___ += fz;
            atom_i.F_nonbonded_z___ -= fz;
        }
    }
}

}//MM

