#ifdef WIN32
#pragma warning(disable : 4786)
#endif

#include "Conventional_interactions.h"

#include "Program_title.h"

#include "Common_interactions.h"
#include "Rough_hydrophobic_interactions.h"

#include "List_interactions.h"
#include "Commonplace_interactions.h"
#include "Yakub_Ronchi_interactions.h"
//#include "Still_interactions.h"

#ifndef ASCALAPH
#include "Solvation_interactions.h"
#endif

//#include "CUDA_interactions.h"

#include "Log.h"

#include "Create.h"
#include "User.h"
#include "Geometry.h"
#include "Defs.h"
#include "MOT.h"

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

#include "Point_3D_impl.h"
#include "Dihedral_angle.h"

#include <math.h>

extern "C" bool GPU_is_on;
bool GPU_is_on = false;

namespace MM
{
const int Conventional_interactions::water = ::mot().atoi ("OW");


//Common_interactions * GPU;
Common_interactions * interactions_prototype = 0;

//==============================================================================
Common_interactions * Common_interactions::
create (Model & model/*, Force_field & ff*/)
{
#ifdef ASCALAPH
    if (interactions_prototype == 0)
    {
        //log() << "\nYakub_Ronchi_interactions\n";
        //return new Yakub_Ronchi_interactions (model/*, ff*/);
        return new List_interactions (model/*, ff*/);
    }
    else
    {
        //log() << "\nMay be CUDA_interactions\n";
        return interactions_prototype->clone (model/*, ff*/);
    }
#else
    //return new Yakub_Ronchi_interactions (model/*, ff*/);
    //log() << "\nSolvation_interactions\n";
    return new Solvation_interactions (model/*, ff*/);
    //return new Still_interactions (model/*, ff*/);

#endif
}

//Common_interactions * Common_interactions::
//create (Atom_group & atom_group, Bond_group & bond_group, Force_field & ff)
//{
////    return new Commonplace_interactions (atom_group, bond_group, ff);
//    //return new Rough_hydrophobic_interactions (atom_group, bond_group, ff);
//
//    //return new Conventional_interactions (atom_group, bond_group, ff);
//    //return new List_interactions (atom_group, bond_group, ff);
//
//    //if (GPU_support)
    //    //    return new CUDA_interactions (atom_group, bond_group, ff);
//    //else
//    //    return new Yakub_Ronchi_interactions (atom_group, bond_group, ff);
//
//    if (interactions_prototype == 0)
//        return new Yakub_Ronchi_interactions (atom_group, bond_group, ff);
//    else
//        return interactions_prototype->clone (atom_group, bond_group, ff);
//}
//==============================================================================


double Conventional_interactions::
potential ()
{
    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.;

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

    try
    {
        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 ("Common_interactions: Unknown exception.");
    }

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

    return last_potential_;
}

void Conventional_interactions::
zero_forces ()
{
    int count = atom_count ();

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

        atom.force_x___ = 0.;
        atom.force_y___ = 0.;
        atom.force_z___ = 0.;

        atom.F_bond_stretch_x___ = 0.;
        atom.F_bond_stretch_y___ = 0.;
        atom.F_bond_stretch_z___ = 0.;

        atom.F_angle_x___ = 0.;
        atom.F_angle_y___ = 0.;
        atom.F_angle_z___ = 0.;

        atom.F_torsion_x___ = 0.;
        atom.F_torsion_y___ = 0.;
        atom.F_torsion_z___ = 0.;

        atom.F_improper_x___ = 0.;
        atom.F_improper_y___ = 0.;
        atom.F_improper_z___ = 0.;

        atom.F_nonbonded_x___ = 0.;
        atom.F_nonbonded_y___ = 0.;
        atom.F_nonbonded_z___ = 0.;

        //atom.FS_nonbonded_x___ = 0.;
        //atom.FS_nonbonded_y___ = 0.;
        //atom.FS_nonbonded_z___ = 0.;

        //atom.FL_nonbonded_x___ = 0.;
        //atom.FL_nonbonded_y___ = 0.;
        //atom.FL_nonbonded_z___ = 0.;
    }

    //set_virial (0,0,0,0);
}

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

    int i, count;

    //profile().interactions___add_force___others ();

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

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


    try
    {
        count = atom_count ();

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

        profile().interactions___add_force___others ();
        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();

        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();
    }
    catch(...)
    {
        to_user().fatal_error ("Common_interactions: Unknown exception.");
    }

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

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

/*void Nonbonded_interactions::
potential (Atom_cache const & atom)
{
    //__asm nop

    double  lIepsilon = 1. / epsilon_;
    double  A = 0., B = 0.;

    int     atom_number = atom.number___;
    int     count       = atom_count();

    double  charge_332  = atom.charge___ * EPS0; //331.8;
    double  e           = 0;
    double  x, y, z;
    double  lIr_6, lIr_eps;

    Force_field::VdWaals * vdw = force_field_.VdW_par [atom.mm_type___];

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

    miss_neighbours (atom);

    for (int i=atom_number+1;  i<count;  ++i)
    {
        Atom_cache & current_atom = atom_[i];

        Force_field::VdWaals & tmp = vdw [current_atom.mm_type___];  //!!!
        A = tmp.A;
        B = tmp.B;    

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

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

        e += current_atom.transparency___ * (
            (A * lIr_6 - B ) * lIr_6 
            + current_atom.charge___ * charge_332 * lIr_eps);
    }

    process_14 (atom, false);
    restore_neighbours ();

    last_potential_ += e;
}//*/

void Conventional_interactions::
add_atom_force (Atom_cache & atom)
{
    //Timer timer;    Text  time;
    //profile().interactions___add_force___others ();
    TIME ("add_atom_force {              ");

    Force_field & ff = force_field();

    add_nonbonded_force        (atom);  
    //time += "   nonbonded_force "; time += timer.elapsed();    timer.restart();
    TIME ("nonbonded_force               ");

    profile().interactions___add_force___add_full_nonbonded_force ();

    add_angle_force            (atom, ff);
    //time += "   angle_force "; time += timer.elapsed();    timer.restart();
    TIME ("angle_force                   ");

    add_improper_torsion_force (atom);
    //time += "   improper_torsion_force "; time += timer.elapsed();    timer.restart();
    TIME ("improper_torsion_force        ");

    //time += " sec\n";
    //log () << time;
    TIME ("add_atom_force }              ");
    //profile().interactions___add_force___add_atom_force ();
}

void Conventional_interactions::
add_nonbonded_force (Atom_cache & atom)
{
    double  lIepsilon = 1. / epsilon();
    double  A = 0., B = 0.;

    int     atom_number = atom.number___;
    int     count       = atom_count();

    double  charge_332  = atom.charge___ * EPS0; //331.8;
    double  x, y, z;
    double  lIr_2, lIr_6, lIr_eps;

    Force_field::VdWaals * vdw = force_field().VdW_par [atom.mm_type___];

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

    miss_neighbours (atom);

    for (int i=atom_number+1;  i<count;  ++i)
    {
        Atom_cache & current_atom = atom_[i];

        Force_field::VdWaals & tmp = vdw [current_atom.mm_type___];
        A = tmp.A;
        B = tmp.B;    

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

        lIr_2 = 1.0 / (x*x + y*y + z*z);

        lIr_eps = sqrt (lIr_2) * lIepsilon;
        //lIr_eps = 2.0 * lIr_2 * lIepsilon;
        
        lIr_6 = lIr_2 * lIr_2 * lIr_2;

        double q = current_atom.charge___ * charge_332 * lIr_eps;
        double K = - current_atom.transparency___ * 
            (6 * (B - 2 * A * lIr_6) * lIr_6 - q) * lIr_2;

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

        atom.        force_x___ += fx;
        current_atom.force_x___ -= fx;

        atom.        force_y___ += fy;
        current_atom.force_y___ -= fy;

        atom.        force_z___ += fz;
        current_atom.force_z___ -= fz;
    }
    process_14 (atom, Force);
    restore_neighbours ();
}//*/

void Conventional_interactions::
add_angle_force (Atom_cache & atom_j, Force_field & ff)
{
    if (atom_j.mm_type___ == water)
    {
        SPCf_potential (atom_j, true);
        return;
    }

    double  K=0,  angleEq=0;
    int     number_of_bond = atom_j.bound_count();
    int     i, k;
    //Force_field & ff = force_field();

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

        for (k=i+1;  k<number_of_bond;  ++k)
        {
            Atom_cache & 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;
            }


            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 Conventional_interactions::
add_improper_torsion_force (Atom_cache & atom_k)
{
    if (atom_k.bound_count() == 3)
    {
        Atom_cache & atom_i = atom_k.bound(0);
        Atom_cache & atom_j = atom_k.bound(1);
        Atom_cache & atom_l = atom_k.bound(2);

        const int mm_type_i = atom_i.mm_type___;
        const int mm_type_j = atom_j.mm_type___;
        const int mm_type_k = atom_k.mm_type___;
        const int mm_type_l = atom_l.mm_type___;

        int i, j, l;
        const Force_field::Imp_Tor *improper_tortion = force_field().imp_tor
            (mm_type_i, mm_type_j, mm_type_k, mm_type_l,   i, j, l);

        if (improper_tortion != 0) //fix ijl order
        {
            const double Phi0 = __Pi__; //__Pi__180.
            const double V_2  = improper_tortion->half_V;
            const double n    = improper_tortion->N;

            double Phi = MM::torsion (atom_i, atom_j, atom_k, atom_l); 
            double dV_I_dPhi = - V_2 * n * sin (n * Phi - Phi0);

            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);

            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_improper_x___ += Fi.x(); 
            atom_i.F_improper_y___ += Fi.y(); 
            atom_i.F_improper_z___ += Fi.z();

            atom_l.F_improper_x___ += Fl.x(); 
            atom_l.F_improper_y___ += Fl.y(); 
            atom_l.F_improper_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_improper_x___ += Fj.x(); 
            atom_j.F_improper_y___ += Fj.y(); 
            atom_j.F_improper_z___ += Fj.z();

            atom_k.F_improper_x___ += Fk.x(); 
            atom_k.F_improper_y___ += Fk.y(); 
            atom_k.F_improper_z___ += Fk.z();
        }
    }
}

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

void Conventional_interactions::
add_bond_stretch_force (Bond_cache & bond)
{
    // Harmonic Bond =========================================================
    double Kr, Req;

    Atom_cache & atom_1 = *bond.first_atom___;
    Atom_cache & atom_2 = *bond.second_atom___;

    if (atom_1.mm_type___ == water || atom_2.mm_type___ == water)
        return;

    force_field().bond (atom_1.mm_type___, atom_2.mm_type___,  Kr, Req );

    if (ad_hoc() && force_field().fail())
    {
        Kr = 300.0;
        if (atom_1.origin___->element().number() == 1 || 
            atom_2.origin___->element().number() == 1)
            Req = 1;
        else
            Req = 1.5;
    }
    double dx = atom_2.x___ - atom_1.x___;
    double dy = atom_2.y___ - atom_1.y___;
    double dz = atom_2.z___ - atom_1.z___;

    double R           = distance (atom_1, atom_2);
    double R_minus_Req = R - Req;
    double dPot_I_dR   = 2. * Kr * R_minus_Req;

    double dPot_I_dR_I_R =  dPot_I_dR / R;
    double fx = dx * dPot_I_dR_I_R;
    double fy = dy * dPot_I_dR_I_R;
    double fz = dz * dPot_I_dR_I_R;

    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_bond_stretch_x___ += fx;
    atom_2.F_bond_stretch_x___ -= fx;

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

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

    add_virial (-dPot_I_dR_I_R*R*R,
                -dPot_I_dR_I_R*dx*dx,
                -dPot_I_dR_I_R*dy*dy,
                -dPot_I_dR_I_R*dz*dz);
}

void Conventional_interactions::
add_torsion_force (Bond_cache & bond_j)
{
    Atom_cache & atom_j = *bond_j.first_atom___;
    Atom_cache & atom_k = *bond_j.second_atom___;

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

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

        //for (j=0;  j<atom_3.bound_count();  j++)
        for (l=0;  l<atom_k.bound_count();  ++l)
        {
            //Atom_cache const & atom_4 = atom_3.bound (j);
            Atom_cache & atom_l = atom_k.bound (l);
            //if (&atom_j == &atom_4) continue;
            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;

            //k_max = tor_tmp->size();
            //for (k=0;  k<k_max;  k++)
            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;

                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();

                ////    Fi = getTorsion( A1, A2, A3, A4 );
                //double  ab_x, ab_y, ab_z;
                //ab_x= atom_2.x___ -  atom_1.x___;
                //ab_y= atom_2.y___ -  atom_1.y___;
                //ab_z= atom_2.z___ -  atom_1.z___;
                //double  bc_x, bc_y, bc_z;
                //bc_x= atom_3.x___ -  atom_j.x___;
                //bc_y= atom_3.y___ -  atom_j.y___;
                //bc_z= atom_3.z___ -  atom_j.z___;
                //double  cd_x, cd_y, cd_z;
                //cd_x= atom_4.x___ -  atom_3.x___;
                //cd_y= atom_4.y___ -  atom_3.y___;
                //cd_z= atom_4.z___ -  atom_3.z___;

                ////    TVector  n1 = (ab*bc)();
                ////    TVector  n2 = (bc*cd)();
                //double  Xn1, Yn1, Zn1,  Xn2, Yn2, Zn2,  Norm;
                //Xn1 = ab_y * bc_z - ab_z * bc_y;
                //Yn1 = ab_z * bc_x - ab_x * bc_z;
                //Zn1 = ab_x * bc_y - ab_y * bc_x;

                //Xn2 = bc_y * cd_z - bc_z * cd_y;
                //Yn2 = bc_z * cd_x - bc_x * cd_z;
                //Zn2 = bc_x * cd_y - bc_y * cd_x;

                //Norm = sqrt ((Xn1*Xn1+Yn1*Yn1+Zn1*Zn1)*
                //             (Xn2*Xn2+Yn2*Yn2+Zn2*Zn2));

                ////    real cos_nFi = n1^n2;
                //double cos_nFi = (Xn1*Xn2+Yn1*Yn2+Zn1*Zn2) / Norm;

                //switch (n)                                    //fix  for Fi0 0 or 180 only
                //{
                //    case 1:   e += (Fi0 == 0) ?
                //                    V_2 * (1 + cos_nFi) / Mult :
                //                    V_2 * (1 - cos_nFi) / Mult;
                //              break;

                //    case 2:   e += (Fi0 == 0) ?
                //                    V_2 * (2*cos_nFi*cos_nFi) / Mult :
                //                    V_2 * 2 * (1 - cos_nFi*cos_nFi) / Mult;
                //              break;

                //    case 3:   e += (Fi0 == 0) ?
                //                    V_2 * (1 + (4*cos_nFi*cos_nFi*cos_nFi -
                //                    3 * cos_nFi)) / Mult :
                //                    V_2 * (1 - (4*cos_nFi*cos_nFi*cos_nFi -
                //                    3 * cos_nFi)) / Mult;
                //              break;

                //    default:  real  Fi = MM::torsion
                //                (atom_1, atom_2, atom_3, atom_4);

                //              e += V_2 * (1 + cos (n*Fi - Fi0)) / Mult;
                //}
            }
        }                
    }
}

void Conventional_interactions::
atom_potential (Atom_cache & atom)
{
    nonbonded_potential        (atom);
    angle_potential            (atom);
    improper_torsion_potential (atom);
}

void Conventional_interactions::
bond_potential (Bond_cache const & bond)
{
    bond_stretch_potential (bond);
    torsion_potential (bond);
}

void Conventional_interactions::
nonbonded_potential (Atom_cache & atom)
{
    // Nonbonded interactions =================================================

    double  lIepsilon = 1. / epsilon();
    double  A = 0., B = 0.;

    int     atom_number = atom.number___;
    int     count       = atom_count();

    double  charge_332  = atom.charge___ * EPS0; //331.8;
    //double  e           = 0;
    double  P_VdW = 0., P_el = 0.;
    double  x, y, z;
    double  lIr_6, lIr_eps;

    Force_field::VdWaals * vdw = force_field().VdW_par [atom.mm_type___];

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

    miss_neighbours (atom);

    // 0.00315
    //for (int i=atom_number+1;  i<count;  i+=4)
    for (int i=atom_number+1;  i<count;  ++i)
    {
        Atom_cache & current_atom  = atom_[i];
        //Atom_cache & current_atom2 = atom_[i+1];
        //Atom_cache & current_atom3 = atom_[i+2];
        //Atom_cache & current_atom4 = atom_[i+3];

        Force_field::VdWaals & tmp = vdw [current_atom.mm_type___];  //!!!
        A = tmp.A;
        B = tmp.B; 

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

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

        //e += current_atom.transparency___ * (
        //    (A * lIr_6 - B ) * lIr_6 
        //    + current_atom.charge___ * charge_332 * lIr_eps);
        P_VdW += current_atom.transparency___ * 
            (A * lIr_6 - B ) * lIr_6;

        P_el += current_atom.transparency___ *
            current_atom.charge___ * charge_332 * lIr_eps;

        //Force_field::VdWaals & tmp2 = vdw [current_atom2.mm_type___];  //!!!
        //A = tmp2.A;
        //B = tmp2.B;    

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

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

        ////e += current_atom2.transparency___ * (
        ////    (A * lIr_6 - B ) * lIr_6 
        ////    + current_atom2.charge___ * charge_332 * lIr_eps);
        //P_VdW += current_atom2.transparency___ * 
        //    (A * lIr_6 - B ) * lIr_6;

        //P_el += current_atom2.transparency___ *
        //    current_atom2.charge___ * charge_332 * lIr_eps;


        //Force_field::VdWaals & tmp3 = vdw [current_atom3.mm_type___];  //!!!
        //A = tmp3.A;
        //B = tmp3.B;    

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

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

        ////e += current_atom3.transparency___ * (
        ////    (A * lIr_6 - B ) * lIr_6 
        ////    + current_atom3.charge___ * charge_332 * lIr_eps);
        //P_VdW += current_atom3.transparency___ * 
        //    (A * lIr_6 - B ) * lIr_6;

        //P_el += current_atom3.transparency___ *
        //    current_atom3.charge___ * charge_332 * lIr_eps;

        //Force_field::VdWaals & tmp4 = vdw [current_atom4.mm_type___];  //!!!
        //A = tmp4.A;
        //B = tmp4.B;    

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

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

        ////e += current_atom4.transparency___ * (
        ////    (A * lIr_6 - B ) * lIr_6 
        ////    + current_atom4.charge___ * charge_332 * lIr_eps);
        //P_VdW += current_atom4.transparency___ * 
        //    (A * lIr_6 - B ) * lIr_6;

        //P_el += current_atom4.transparency___ *
        //    current_atom4.charge___ * charge_332 * lIr_eps;
    }

    process_14 (atom, Potential);
    restore_neighbours ();

    Van_der_Waals_potential_ += P_VdW;
    electrostatic_potential_ += P_el;
    last_potential_ += P_el + P_VdW;
}

void Conventional_interactions::
angle_potential (Atom_cache const & atom_2)
{
    if (atom_2.mm_type___ == water)
    {
        SPCf_potential (const_cast<Atom_cache &>(atom_2));
        return;
    }

    double  e = 0,  K=0,  angleEq=0;
    Force_field & ff = force_field();

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

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

        for (k=i+1;  k<number_of_bond;  k++)
        {
            Atom_cache const & atom_3 = atom_2.bound (k);
            ff.angle (atom_1.mm_type___, atom_2.mm_type___, atom_3.mm_type___,
                K, angleEq );

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

            double  angle_minus_angleEq = 
                MM::angle (atom_1, atom_2, atom_3) - angleEq;

            e += K * angle_minus_angleEq * angle_minus_angleEq;
        }
    }
    angle_potential_ += e;
    last_potential_  += e; 
}

void Conventional_interactions::
SPCf_potential (Atom_cache & atom_2, bool force)
{
    //if (force)
    //    profile().interactions___add_force___others ();

    double OH_limit = 1.3; // A
    //double OH_limit = 1000; // A
    //double HH_limit = 2.0; // A

    // dyn  = g sm / sec2
    // N    = kg m / sec2
    // dyn  = 1e-5 N
    // mdyn = 1e-8 N
    
    // AKMA, kcal/mol-A = 6.95 10-1 N = 6.95 x 10+4 dyn = 6.95 x 10+7 mdyn
    
    // kJ/mol-A = 8.987551787 e16 / 6.0221367 e23 =  54.124265
    // kg = 2.0614841 e34 Hartree * 627.51 = 1293.601887591 e34 kcal mol-1

    // factor 602.19912

    double  u = Constant::kJ_to_kcal; 

    //                kJ/mol-A      mdyn/A
    double  a   = u * 5619.12;  //  9.331;       // 9.331 mdyn/A = 9.331 / 6.95 x 10+7 = 1.342e-7    = 5.618e-7 kJ/mol-A 
    //double  aI2 = u * 2809.56;  //  9.331 / 2.;                                                      = 2.8088e-7 kJ/mol-A 
    double  b   = u * 1374.82;  //  2.283;
    double  bI2 = u *  687.41;  //  2.283;
    double  c   = u * -884.63;  // -1.469; 
    double  d   = u *  467.31;  //  0.776; 
    //double  D = u * 426.35698; //  0.708;
    
    double  p = 2.566; // 1/A
    // 2 p^2 D = a      5619.12 / 2 * 2.566^2 = 426.70 kJ/mol-A 
    double D_OH = a / (2 * p*p);

    double K_quadratic = 
        ((1 - exp(-p * (OH_limit-1)))*(1 - exp(-p * (OH_limit-1)))) 
        / ((OH_limit-1) * (OH_limit-1));

    Atom_cache & atom_1 = atom_2.bound (0);
    Atom_cache & atom_3 = atom_2.bound (1);

    double R1_OH_x = atom_1.x___ - atom_2.x___;
    double R1_OH_y = atom_1.y___ - atom_2.y___;
    double R1_OH_z = atom_1.z___ - atom_2.z___;

    double R2_OH_x = atom_3.x___ - atom_2.x___;
    double R2_OH_y = atom_3.y___ - atom_2.y___;
    double R2_OH_z = atom_3.z___ - atom_2.z___;

    double R3_HH_x = atom_3.x___ - atom_1.x___;
    double R3_HH_y = atom_3.y___ - atom_1.y___;
    double R3_HH_z = atom_3.z___ - atom_1.z___;


    double R1_OH = sqrt (R1_OH_x*R1_OH_x + R1_OH_y*R1_OH_y + R1_OH_z*R1_OH_z);
    double R2_OH = sqrt (R2_OH_x*R2_OH_x + R2_OH_y*R2_OH_y + R2_OH_z*R2_OH_z);
    double R3_HH = sqrt (R3_HH_x*R3_HH_x + R3_HH_y*R3_HH_y + R3_HH_z*R3_HH_z);

    //fix
    //R1_OH = R1_OH < OH_limit ? R1_OH : OH_limit;
    //R2_OH = R2_OH < OH_limit ? R2_OH : OH_limit;
    //R3_HH = R3_HH < HH_limit ? R3_HH : HH_limit;

    double F1, F2, E1, E2;

    double dR1_OH = R1_OH - 1.;           // Reqv = 1 A
    double dR2_OH = R2_OH - 1.;
    //double dR3_HH = R3_HH - 1.63328311;   // Reqv = 1.633 A
    double dR3_HH = R3_HH - 1.633;   // Reqv = 1.63328311 A
    
    double exp_1 = exp (-p * dR1_OH);
    double exp_2 = exp (-p * dR2_OH);
    
    double I_exp_1 = 1 - exp_1;
    double I_exp_2 = 1 - exp_2;

    if (force)
    {
        if (R1_OH < OH_limit)
            F1 = - 2. * D_OH * I_exp_1 * exp_1 * p - c * dR3_HH - d * dR2_OH;
        else
            F1 = - K_quadratic * 2. * dR1_OH * D_OH * p - c * dR3_HH - d * dR1_OH;

        if (R2_OH < OH_limit)
            F2 = - 2. * D_OH * I_exp_2 * exp_2 * p - c * dR3_HH - d * dR1_OH;
        else
            F2 = - K_quadratic * 2. * dR2_OH * D_OH * p - c * dR3_HH - d * dR1_OH;

        double F3 = - b * dR3_HH - c * (dR1_OH + dR2_OH); 

        F1 /= R1_OH;    F2 /= R2_OH;   F3 /= R3_HH;

        double f1x = F1 * R1_OH_x;
        double f1y = F1 * R1_OH_y;
        double f1z = F1 * R1_OH_z;

        double f2x = F2 * R2_OH_x;
        double f2y = F2 * R2_OH_y;
        double f2z = F2 * R2_OH_z;

        double f3x = F3 * R3_HH_x;
        double f3y = F3 * R3_HH_y;
        double f3z = F3 * R3_HH_z;

        atom_2.force_x___ += -f1x - f2x;
        atom_2.force_y___ += -f1y - f2y;
        atom_2.force_z___ += -f1z - f2z;

        atom_1.force_x___ +=  f1x - f3x;
        atom_1.force_y___ +=  f1y - f3y;
        atom_1.force_z___ +=  f1z - f3z;

        atom_3.force_x___ +=  f2x + f3x;
        atom_3.force_y___ +=  f2y + f3y;
        atom_3.force_z___ +=  f2z + f3z;

        atom_2.F_bond_stretch_x___ += -f1x - f2x;
        atom_2.F_bond_stretch_y___ += -f1y - f2y;
        atom_2.F_bond_stretch_z___ += -f1z - f2z;

        atom_1.F_bond_stretch_x___ +=  f1x - f3x;
        atom_1.F_bond_stretch_y___ +=  f1y - f3y;
        atom_1.F_bond_stretch_z___ +=  f1z - f3z;

        atom_3.F_bond_stretch_x___ +=  f2x + f3x;
        atom_3.F_bond_stretch_y___ +=  f2y + f3y;
        atom_3.F_bond_stretch_z___ +=  f2z + f3z;
    
        add_virial ((F1* R1_OH*R1_OH + F2* R2_OH*R2_OH + F3* R3_HH*R3_HH), 
                    0, 0, 0);
    }
    else
    {
        if (R1_OH < OH_limit)
            E1 = D_OH * I_exp_1*I_exp_1 + c * dR1_OH * dR3_HH;
        else
            E1 = K_quadratic * D_OH * dR1_OH*dR1_OH  + c * dR1_OH * dR3_HH;

        if (R2_OH < OH_limit)
            E2 = D_OH * I_exp_2*I_exp_2 + c * dR2_OH * dR3_HH;
        else
            E2 = K_quadratic * D_OH * dR2_OH*dR2_OH  + c * dR2_OH * dR3_HH;

        double E3 = bI2 * dR3_HH * dR3_HH  + d * dR1_OH * dR2_OH;

        double e = E1 + E2 + E3;
        water_potential_ += e;
        last_potential_  += e; 
    }
    //if (force)
    //    profile().interactions___add_force___add_water_force ();
}

void Conventional_interactions::
improper_torsion_potential  (Atom_cache const & atom_k)
{
    if (atom_k.bound_count() == 3)
    {
        Atom_cache const & atom_i = atom_k.bound(0);
        Atom_cache const & atom_j = atom_k.bound(1);
        Atom_cache const & atom_l = atom_k.bound(2);

        const int mm_type_i = atom_i.mm_type___;
        const int mm_type_j = atom_j.mm_type___;
        const int mm_type_k = atom_k.mm_type___;
        const int mm_type_l = atom_l.mm_type___;

        int i, j, l;
        const Force_field::Imp_Tor *improper_tortion = force_field().imp_tor
            (mm_type_i, mm_type_j, mm_type_k, mm_type_l,   i, j, l);

        if (improper_tortion != 0) //fix ijl order
        {
            const double Phi0 = __Pi__; //__Pi__180.
            const double V_2  = improper_tortion->half_V;
            const double n    = improper_tortion->N;

            double Phi = MM::torsion (atom_i, atom_j, atom_k, atom_l); 
            double e = V_2 * (1 + cos (n*Phi - Phi0));
            
            improper_torsion_potential_ += e;
            last_potential_ += e;
        }
    }
}

void Conventional_interactions::
bond_stretch_potential (Bond_cache const & bond)
{
    // Harmonic Bond =========================================================
    double Kr, Req;

    Atom_cache const & atom_1 = *bond.first_atom___;
    Atom_cache const & atom_2 = *bond.second_atom___;

    if (atom_1.mm_type___ == water || atom_2.mm_type___ == water)
        return;

    force_field().bond (atom_1.mm_type___, atom_2.mm_type___,  Kr, Req );

    if (ad_hoc() && force_field().fail())
    {
        Kr = 300.0;
        if (atom_1.origin___->element().number() == 1 || 
            atom_2.origin___->element().number() == 1)
            Req = 1;
        else
            Req = 1.5;
    }

    double  R_minus_Req = distance (atom_1, atom_2) - Req;

    double e = Kr * R_minus_Req * R_minus_Req;
    
    bond_potential_ += e;
    last_potential_ += e;
}

void Conventional_interactions::
torsion_potential (Bond_cache const & bond)
{
    // e = V/2 ( 1 + cos( n*Fi - Fi0 ) )
    double e = 0;
    const Atom_cache & atom_2 = *bond.first_atom___;
    const Atom_cache & atom_3 = *bond.second_atom___;

    int     i,j,k, k_max;
    Force_field & ff = force_field();

    for (i=0;  i<atom_2.bound_count();  i++)
    {
        Atom_cache const & atom_1 = atom_2.bound (i);
        if (&atom_1 == &atom_3) continue;

        for (j=0;  j<atom_3.bound_count();  j++)
        {
            Atom_cache const & atom_4 = atom_3.bound (j);
            if (&atom_2 == &atom_4) continue;
                                                           //fix may be 0             
            const Array <Force_field::Torsions> * tor_tmp; //!!!
            tor_tmp = ff.torsion (atom_1.mm_type___, 
                                            atom_2.mm_type___, 
                                            atom_3.mm_type___, 
                                            atom_4.mm_type___);

            //fix
            if (tor_tmp == 0)
                 continue;

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

                if (V_2 < d_mach_eps_2)
                    continue;

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

                //    Fi = getTorsion( A1, A2, A3, A4 );
                double  ab_x, ab_y, ab_z;
                ab_x= atom_2.x___ -  atom_1.x___;
                ab_y= atom_2.y___ -  atom_1.y___;
                ab_z= atom_2.z___ -  atom_1.z___;
                double  bc_x, bc_y, bc_z;
                bc_x= atom_3.x___ -  atom_2.x___;
                bc_y= atom_3.y___ -  atom_2.y___;
                bc_z= atom_3.z___ -  atom_2.z___;
                double  cd_x, cd_y, cd_z;
                cd_x= atom_4.x___ -  atom_3.x___;
                cd_y= atom_4.y___ -  atom_3.y___;
                cd_z= atom_4.z___ -  atom_3.z___;

                //    TVector  n1 = (ab*bc)();
                //    TVector  n2 = (bc*cd)();
                double  Xn1, Yn1, Zn1,  Xn2, Yn2, Zn2,  Norm;
                Xn1 = ab_y * bc_z - ab_z * bc_y;
                Yn1 = ab_z * bc_x - ab_x * bc_z;
                Zn1 = ab_x * bc_y - ab_y * bc_x;

                Xn2 = bc_y * cd_z - bc_z * cd_y;
                Yn2 = bc_z * cd_x - bc_x * cd_z;
                Zn2 = bc_x * cd_y - bc_y * cd_x;

                Norm = sqrt ((Xn1*Xn1+Yn1*Yn1+Zn1*Zn1)*
                             (Xn2*Xn2+Yn2*Yn2+Zn2*Zn2)) + d_mach_eps_2;

                //    real cos_nFi = n1^n2;
                double cos_nFi = (Xn1*Xn2+Yn1*Yn2+Zn1*Zn2) / Norm;

                switch (n)                                    //fix  for Fi0 0 or 180 only
                {
                    case 1:   e += (Fi0 == 0) ?
                                    V_2 * (1 + cos_nFi) / Mult :
                                    V_2 * (1 - cos_nFi) / Mult;
                              break;

                    case 2:   e += (Fi0 == 0) ?
                                    V_2 * (2*cos_nFi*cos_nFi) / Mult :
                                    V_2 * 2 * (1 - cos_nFi*cos_nFi) / Mult;
                              break;

                    case 3:   e += (Fi0 == 0) ?
                                    V_2 * (1 + (4*cos_nFi*cos_nFi*cos_nFi -
                                    3 * cos_nFi)) / Mult :
                                    V_2 * (1 - (4*cos_nFi*cos_nFi*cos_nFi -
                                    3 * cos_nFi)) / Mult;
                              break;

                    default:  real  Fi = MM::torsion
                                (atom_1, atom_2, atom_3, atom_4);

                              e += V_2 * (1 + cos (n*Fi - Fi0)) / Mult;
                }
            }
        }                
    }

    torsion_potential_ += e;
    last_potential_    += e;
}

/*void Nonbonded_interactions::
add_force (Atom_cache & atom)
{
    double  lIepsilon = 1. / epsilon_;
    double  A = 0., B = 0.;

    int     atom_number = atom.number___;
    int     count       = atom_count();

    double  charge_332  = atom.charge___ * EPS0; //331.8;
    double  x, y, z;
    double  lIr_2, lIr_6, lIr_eps;

    Force_field::VdWaals * vdw = force_field_.VdW_par [atom.mm_type___];

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

    miss_neighbours (atom);

    for (int i=atom_number+1;  i<count;  i+=2)
    {
        Atom_cache & current_atom  = atom_[i];
        Atom_cache & current_atom2 = atom_[i+1];

        Force_field::VdWaals & tmp = vdw [current_atom.mm_type___];
        A = tmp.A;
        B = tmp.B;    

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

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

        double q = current_atom.charge___ * charge_332 * lIr_eps;
        double K = - current_atom.transparency___ * 
            (6 * (B - 2 * A * lIr_6) * lIr_6 - q) * lIr_2;

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

        atom.        force_x___ += fx;
        current_atom.force_x___ -= fx;

        atom.        force_y___ += fy;
        current_atom.force_y___ -= fy;

        atom.        force_z___ += fz;
        current_atom.force_z___ -= fz;


        Force_field::VdWaals & tmp2 = vdw [current_atom2.mm_type___];
        A = tmp2.A;
        B = tmp2.B;    

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

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

        q = current_atom2.charge___ * charge_332 * lIr_eps;
        K = - current_atom2.transparency___ * 
            (6 * (B - 2 * A * lIr_6) * lIr_6 - q) * lIr_2;

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

        atom.         force_x___ += fx;
        current_atom2.force_x___ -= fx;

        atom.         force_y___ += fy;
        current_atom2.force_y___ -= fy;

        atom.         force_z___ += fz;
        current_atom2.force_z___ -= fz;
    }

    restore_neighbours ();
}//*/

void Conventional_interactions::
process_14 (Atom_cache & atom_1, Task task)
{
    for (int i=0;  i<atom_1.bound_count();  i++)
    {
        Atom_cache & atom_2 = atom_1.bound (i);

        for (int j=0;  j<atom_2.bound_count();  j++)
        {
            Atom_cache & atom_3 = atom_2.bound (j);
            if (&atom_1 == &atom_3) continue;

            for (int k=0;  k<atom_3.bound_count();  k++)
            {
                Atom_cache & atom_4 = atom_3.bound (k);
                
                if (&atom_2 == &atom_4 || 
                     atom_4.flag___.is (Atom_kit::missed)) 
                        continue;
                
                atom_4.flag___.set (Atom_kit::missed);
                miss_set_.insert (&atom_4);

                if (atom_1.number___ >= atom_4.number___)
                    continue;

                switch (task)
                {
                case Force:     add_nonbonded_force_14 (atom_1, atom_4); break;
                case Potential: nonbonded_potential_14 (atom_1, atom_4); break;
                case List:      list_nonbonded_14      (atom_1, atom_4);  break;
                case No:        break;
                default:        FLAW ("Unknown task.");
                }
            }
        }
    }
}

void Conventional_interactions::
nonbonded_potential_14 (Atom_cache & atom_1, Atom_cache & atom_2)
{
//    double  e = 0;
    double  P_VdW = 0., P_el = 0.;
    double  A, B;
    Force_field & ff = force_field();

    // e = A/r^12 - B/r^6  +  electrostatics
    ff.VdW (atom_1.mm_type___, atom_2.mm_type___, A, 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 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

    P_VdW += (A * lIr_6 - B) * lIr_6 * ff.VdW_Scale;
    P_el += atom_1.charge___ * atom_2.charge___ * EPS0 * 
        ff.E_Static_Scale * lIr_eps;

    //fix I
    //double R_star_1 = ff.R_eps[atom_1.mm_type___].R_star;
    //double IR_star3_1 = ff.I_ * R_star_1 * R_star_1 * R_star_1;
    //double q2_1 = atom_1.charge___ * atom_1.charge___;
    //double R_star_2 = ff.R_eps[atom_2.mm_type___].R_star;
    //double IR_star3_2 = ff.I_ * R_star_2 * R_star_2 * R_star_2;
    //double q2_2 = atom_2.charge___ * atom_2.charge___;
    //double Ind = -IR_star3_1 * q2_2 * lIr_6 - IR_star3_2 * q2_1 * lIr_6;

    Van_der_Waals_potential_ += P_VdW;
    electrostatic_potential_ += P_el;
    //induction_potential_     += Ind;
    last_potential_          += P_VdW + P_el /*+ Ind*/;
}

void Conventional_interactions::
add_nonbonded_force_14 (Atom_cache & atom_1, Atom_cache & atom_2)
{
    double  A, B;
    Force_field & ff = force_field();

    // e = A/r^12 - B/r^6  +  electrostatics
    ff.VdW (atom_1.mm_type___, atom_2.mm_type___, A, 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)
                  );

    //fix I
    //double R_star_1 = ff.R_eps[atom_1.mm_type___].R_star;
    //double IR_star3_1 = ff.I_ * R_star_1 * R_star_1 * R_star_1;
    //double q2_1 = atom_1.induction_charge___ * atom_1.induction_charge___;
    //double R_star_2 = ff.R_eps[atom_2.mm_type___].R_star;
    //double IR_star3_2 = ff.I_ * R_star_2 * R_star_2 * R_star_2;
    //double q2_2 = atom_2.induction_charge___ * atom_2.induction_charge___;
    //double Ind = -6 * IR_star3_1 * q2_2 * lIr_6 * lIr_2
    //             -6 * IR_star3_2 * q2_1 * lIr_6 * lIr_2;   

    //K += Ind;

    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 Conventional_interactions::
list_nonbonded_14 (Atom_cache &, Atom_cache &)
{
    //fix
}

void Conventional_interactions::
miss_neighbours (Atom_cache const & atom_1)
{
    REQUIRE ("", miss_set_.size() == 0);
    int i, j, k;

    for (i=0;  i<atom_1.bound_count();  i++)
    {
        Atom_cache const & atom_2 = atom_1.bound (i);

        if (atom_2.flag___.is_not (Atom_kit::missed))
        {
            atom_2.transparency___ = 0.;
            atom_2.flag___.set (Atom_kit::missed);
            miss_set_.insert (&atom_2);
        }

        for (j=0;  j<atom_2.bound_count();  j++)
        {
            Atom_cache const & atom_3 = atom_2.bound (j);
            if (&atom_1 == &atom_3) continue;

            if (atom_3.flag___.is_not (MM::Atom_kit::missed))
            {
                atom_3.transparency___ = 0.;
                atom_3.flag___.set (Atom_kit::missed);
                miss_set_.insert (&atom_3);
            }
            
            for (k=0;  k<atom_3.bound_count();  k++)
            {
                Atom_cache const & atom_4 = atom_3.bound (k);

                if (&atom_2 == &atom_4 || &atom_1 == &atom_4) continue;

                if (atom_4.flag___.is_not (Atom_kit::missed))
                    atom_4.transparency___ = 0.;
            }
        }
    }
}

void Conventional_interactions::
restore_neighbours ()
{
    std::set <Atom_cache const *>::iterator it = miss_set_.begin();
    for (;  it!=miss_set_.end();  ++it)
    {
        (*it)->flag___.clear (Atom_kit::missed);
        (*it)->transparency___ = 1.;
    }
    miss_set_.clear();
}

//Conventional_interactions::
//Conventional_interactions (Atom_group & atom_group, Force_field & ff)
//:
//    Common_interactions         (ff),
//    Bond_group_cache (nil<Bond_group>(), atom_group)
//{
//}
//
//Conventional_interactions::
//Conventional_interactions (Atom_group & atom_group, Bond_group & bond_group,
//    Force_field & ff)
//:
//    Common_interactions         (ff),
//    Bond_group_cache (bond_group, atom_group)
//{
//}

Conventional_interactions::
Conventional_interactions (Model & model/*, Force_field & ff*/)
:
    Common_interactions (model/*, ff*/),
    Bond_group_cache    (model, model)
{
}

}//MM

