#include "Yakub_Ronchi_interactions.h"

#include "Model_kit.h"
#include "anint.h"

#include <math.h>

namespace MM
{

//Yakub_Ronchi_interactions::
//Yakub_Ronchi_interactions (Atom_group & atom_group, Force_field & ff)
//:
//    Conventional_interactions (atom_group, ff)
//{
//}
//
//Yakub_Ronchi_interactions::
//Yakub_Ronchi_interactions (Atom_group & atom_group,
//    Bond_group & bond_group, Force_field & ff)
//:
//    Conventional_interactions (atom_group, bond_group, ff)
//{
//}

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

void Yakub_Ronchi_interactions::
nonbonded_potential (Atom_cache & atom_j)
{
    Force_field & ff = force_field();

    double Rm      = rm ();
    double cut_off = cutoff ();
    double I_cut_2 = ff.I_cut_* ff.I_cut_;
    
    double  L = 0., lIL = 0., LI2 = 0.;
//    L = L(box);

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

        if (cut_off > LI2)
        {
            cut_off = LI2;
            if (cutoff_eq_Rm())
                Rm = LI2;
        }
    }

    double cutoff_2 = cut_off*cut_off;

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

    //double  lIepsilon = 1. / epsilon_;//fix
    double  lIepsilon = 1. / force_field().eps;//fix

    double  A = 0., B = 0.;

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

    double  charge_332_j  = atom_j.charge___ * EPS0; //331.8;
    //double  e           = 0;
    double  P_VdW = 0., P_el = 0., P_hph = 0., Yakub_Ronchi = 0.;
    double  x, y, z;
    double  r_2, lIr, lIr_6, lIr_eps;
    
    double  Ind = 0.;                                                           //fix I
    double R_star = ff.R_eps[atom_j.mm_type___].R_star;
    double IR_star3_j = ff.I_ * R_star * R_star * R_star;
    //double q2 = atom.charge___ * atom.charge___;
    double q2_j = atom_j.induction_charge___ * atom_j.induction_charge___;

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

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

    miss_neighbours (atom_j);

    // 0.00315
    for (int i=atom_number+1;  i<count;  ++i)
    //for (int i=atom_number+1;  i<count;  i+=4)
    {
        Atom_cache & atom_i  = atom_[i];

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

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

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

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

        if (r_2 > cutoff_2)
            continue;

        lIr     = 1.0 / sqrt (r_2);
        lIr_eps = lIr * lIepsilon;
        lIr_6   = lIr * lIr;    lIr_6 *= lIr_6 * lIr_6;

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

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

        double current_P_el = atom_i.transparency___ * (r_2 < Rm_2) *
            atom_i.charge___ * charge_332_j * lIr_eps * Yakub_Ronchi;
        P_el += current_P_el;

        //fix I
        //double current_Ind = current_atom.transparency___ * 
        //    -IR_star3 * current_atom.charge___ * current_atom.charge___ * lIr_6;
        double q2_i = atom_i.induction_charge___ * atom_i.induction_charge___;
        double current_Ind = atom_i.transparency___ * 
            -IR_star3_j * (q2_j < I_cut_2) * q2_i * lIr_6;

        double R_star_i = ff.R_eps[atom_i.mm_type___].R_star;
        double IR_star3_i = ff.I_ * R_star_i * R_star_i * R_star_i;
               current_Ind += atom_i.transparency___ * 
            -IR_star3_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;
        }
    }

    process_14 (atom_j, Potential);
    restore_neighbours ();

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

    last_potential_ += P_el + P_VdW + P_hph + Ind;
}

void Yakub_Ronchi_interactions::
add_nonbonded_force (Atom_cache & atom_j)
{
    Force_field & ff = force_field();

    double Rm      = rm ();
    double cut_off = cutoff ();
    double I_cut_2 = ff.I_cut_* ff.I_cut_;

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

        if (cut_off > LI2)
        {
            cut_off = LI2;
            if (cutoff_eq_Rm())
                Rm = LI2;
        }
    }
    
    double cutoff_2 = cut_off*cut_off;

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

    //double  lIepsilon = 1. / epsilon_;
    double  lIepsilon = 1. / force_field().eps;//fix
    double  A = 0., B = 0.;

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

    double  charge_332  = atom_j.charge___ * EPS0; //331.8;   
    //double  charge_332  = atom.charge___ * EPS0 / 2.; //331.8;   

    double  x, y, z;
    double  r_2, lIr, lIr_2, lIr_6/*, lIr_eps*/;

    double R_star = ff.R_eps[atom_j.mm_type___].R_star;
    double IR_star3_j = ff.I_ * R_star * R_star * R_star;
    //double q2 = atom.charge___ * atom.charge___;
    double q2_j = atom_j.induction_charge___ * atom_j.induction_charge___;

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

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

    miss_neighbours (atom_j);

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

        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;

        if (r_2 > cutoff_2)
            continue;

        lIr     = 1.0 / sqrt (r_2);
        lIr_2   = lIr * lIr;
        //lIr_eps = lIr * lIepsilon;
        lIr_6 = lIr_2 * lIr_2 * lIr_2;

        double Kf_Yakub_Ronchi      = lIRm_3     - lIr_2 * lIr; 

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

        double F_el = - lIepsilon * current_atom.charge___ * charge_332 
            * Kf_Yakub_Ronchi * (r_2 < Rm_2);

        //fix I
        //double Ind = -6 * IR_star3 * current_atom.charge___ * current_atom.charge___ * lIr_6 * lIr_2;
        double q2_i = current_atom.induction_charge___ * current_atom.induction_charge___;
        double Ind = -6 * IR_star3_j * (q2_j   < I_cut_2) * q2_i * lIr_6 * lIr_2;
        
        double R_star_i = ff.R_eps[current_atom.mm_type___].R_star;
        double IR_star3_i = ff.I_ * R_star_i * R_star_i * R_star_i;
               Ind += -6 * IR_star3_i * (q2_i < I_cut_2) * q2_j * lIr_6 * lIr_2;

        double K = current_atom.transparency___ * (F_el + F_VdW + Ind);//fix P

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

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

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

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

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

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

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

        add_virial (K*r_2, K*x*x, K*y*y, K*z*z);
    }
    process_14 (atom_j, Force);
    restore_neighbours ();
}

}//MM
