#include "Rough_hydrophobic_interactions.h"

#include <math.h>

namespace MM
{

Rough_hydrophobic_interactions::
Rough_hydrophobic_interactions (Model & model)
:
    Conventional_interactions (model)
{
}

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

//namespace
//{
//    //const double Rm      = 10.;
////    const double Rm      = 5.5;
//    double Rm      = 10.;
//
//    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;
//
//
//
////    const double l       = 0.035;
//    double l       = 0.06;
//    //const double l       = 0.07;
//    //const double l       = 0.09;
//
//    //const double Rm_hph      = 20.;
//    //const double Rm_hph      = 15.;
//    //const double Rm_hph      = 10.;
//    double Rm_hph      = 10.0;
//
//    const double Rm_hph_2    = Rm_hph * Rm_hph;
//    const double lIRm_hph    = 1. / Rm_hph;
//    const double lIRm_hph_3  = 1. / (Rm_hph*Rm_hph*Rm_hph);
//    const double lI2Rm_hph_3 = 0.5 * lIRm_hph_3;
//
//    //const double K_rough = -12.;
//    //const double K_rough = -30;
//    //const double K_rough = -50;
//    //const double K_rough = -1.0;
//    //const double K_rough = -1.;
//    double K_rough = 10;
//}

void Rough_hydrophobic_interactions::
nonbonded_potential (Atom_cache & atom)
{
    double Rm      = force_field().Rm_;
    double l       = force_field().l_;
    double Rm_hph  = force_field().Rm_hph_;
    double K_rough = force_field().K_rough_;

    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;

    const double Rm_hph_2    = Rm_hph * Rm_hph;
    const double lIRm_hph    = 1. / Rm_hph;
    //const double lIRm_hph_3  = 1. / (Rm_hph*Rm_hph*Rm_hph);
    //const double lI2Rm_hph_3 = 0.5 * lIRm_hph_3;

    //const double l      = 0.17;//0.17 -0.774
    //const double l      = 0.06;//0.06 -0.864
    //const double l      = 0.015;//0.015 -0.878
    //const double l      = 0.005;//0.005 -0.826
    //const double l      = 1.25;//0.005 -0.826
    //const double l      = 0.11;//
    //const double l      = 0.035;//

    //const double k      = 1. * (fabs (atom.charge___) - l);
    //const double k      = 1. * (pow(atom.charge___, 2) - l);
    //const double k      = 1. * (pow(atom.charge___, 4) - l);
    //const double k      = 1. * (pow(atom.charge___, 8) - l);
    //const double k      = 1. * (exp(fabs(atom.charge___)) - l);
    //const double k      = 1. * ((pow(atom.charge___, 2)/(pow(atom.charge___, 2)+0.6)) - l);
    const double k      = K_rough * ((pow(atom.charge___, 4)/(pow(atom.charge___, 4)+0.50)) - l);

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

    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., P_hph = 0., Yakub_Ronchi = 0., Yakub_Ronchi_hph = 0;
    double  x, y, z;
    double  r_2, lIr, 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)
    //for (int i=atom_number+1;  i<count;  i+=4)
    {
        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___;

        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;

        //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 rIRm_hph  = r_2 * lIr * lIRm_hph;
        Yakub_Ronchi_hph = 1. + .5 * rIRm_hph * (rIRm_hph * rIRm_hph - 3.); 

        P_VdW += current_atom.transparency___ * 
            (A * lIr_6 - B ) * lIr_6;

        P_el += current_atom.transparency___ * (r_2 < Rm_2) *
            current_atom.charge___ * charge_332 * lIr_eps * Yakub_Ronchi;

        P_hph += current_atom.transparency___ * (r_2 < Rm_hph_2) *
            k * lIr * Yakub_Ronchi_hph;
        //P_hph += current_atom.transparency___ * (r_2 < Rm_2) *
        //    k * lIr * Yakub_Ronchi;
        //P_hph += (1. * (fabs (current_atom.charge___) - l)) * lIr * Yakub_Ronchi;
        //P_hph += (1. * (pow(current_atom.charge___, 2) - l)) * lIr * Yakub_Ronchi;
        //P_hph += (1. * (pow(current_atom.charge___, 4) - l)) * lIr * Yakub_Ronchi;
        //P_hph += (1. * (pow(current_atom.charge___, 8) - l)) * lIr * Yakub_Ronchi;
        //P_hph += (1. * (exp(fabs(current_atom.charge___)) - l)) * lIr * Yakub_Ronchi;
        //P_hph += (1. * ((pow(current_atom.charge___, 2)/(pow(current_atom.charge___, 2)+0.6)) - l)) * lIr * Yakub_Ronchi;
        P_hph += current_atom.transparency___ * (r_2 < Rm_hph_2) *
            (K_rough * ((pow(current_atom.charge___, 4)/(pow(current_atom.charge___, 4) + 0.50)) - l)) * lIr * Yakub_Ronchi_hph;
        //P_hph += current_atom.transparency___ * (r_2 < Rm_2) *
        //    (K_rough * ((pow(current_atom.charge___, 4)/(pow(current_atom.charge___, 4) + 0.50)) - l)) * lIr * Yakub_Ronchi;
    }

    process_14 (atom, Potential);
    restore_neighbours ();

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

    last_potential_ += P_el + P_VdW + P_hph;
}

void Rough_hydrophobic_interactions::
add_nonbonded_force (Atom_cache & atom)
{
    double Rm      = force_field().Rm_;
    double l       = force_field().l_;
    double Rm_hph  = force_field().Rm_hph_;
    double K_rough = force_field().K_rough_;

    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;

    const double Rm_hph_2    = Rm_hph * Rm_hph;
    //const double lIRm_hph    = 1. / Rm_hph;
    const double lIRm_hph_3  = 1. / (Rm_hph*Rm_hph*Rm_hph);
    //const double lI2Rm_hph_3 = 0.5 * lIRm_hph_3;


    const double atom_charge_n = pow(atom.charge___, 4);
    const double k = K_rough * ((atom_charge_n / (atom_charge_n + 0.50)) - l);

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

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

    double  charge_332  = atom.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*/;

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

        r_2     = x*x + y*y + z*z;
        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 Kf_Yakub_Ronchi_hph  = lIRm_hph_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);

        double charge_n = pow (current_atom.charge___, 4);
        double current_k = K_rough * ((charge_n / (charge_n + 0.5)) - l);
        double F_hph = - (k + current_k) * Kf_Yakub_Ronchi_hph * (r_2 < Rm_hph_2);

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

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

}//MM
