#include "Commonplace_interactions.h"

#include "Atom.h"
#include "Atom_group.h"
#include "Atom_kit.h"
#include "Bond.h"
#include "Bond_group.h"
#include "Vector_3D_impl.h"
#include "Atom_numerator.h"
#include "Create.h"

#include <math.h>

namespace MM
{
//Commonplace_interactions::
//Commonplace_interactions (Atom_group & atom_group, Bond_group & bond_group,
//    Force_field & ff)
//:
//    Common_interactions (ff),
//    atom_group_         (atom_group),
//    bond_group_         (bond_group)
//{
//}
//
//Commonplace_interactions::
//Commonplace_interactions (Atom_group & atom_group, 
//    Force_field & ff)
//:
//    Common_interactions (ff),
//    atom_group_         (atom_group),
//    bond_group_         (nil<Bond_group>())
//{
//}

Commonplace_interactions::
Commonplace_interactions (Model & model/*, Force_field & ff*/)
:
    Common_interactions (model/*, ff*/),
    atom_group_         (model),
    bond_group_         (model)
{
}


double Commonplace_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.;

    Atom_numerator numerator(atom_group_);

    try
    {
        count = atom_group_.atom_count();
        for (i=0;  i<count;  ++i)
        {
            Atom & current_atom = atom_group_.atom(i);
            atom_potential (current_atom);
        }

        if (&bond_group_ != &nil<Bond_group>())
        {
            count = bond_group_.bond_count();
            for (i=0;  i<count;  ++i)
            {
                Bond & current_bond = bond_group_.bond(i);
                bond_potential (current_bond);
            }
        }
    }
    catch(...)
    {
        to_user().fatal_error ("Commonplace_interactions: Unknown exception.");
    }

    return last_potential_;
}

void Commonplace_interactions::
zero_forces ()
{
    int             count = atom_group_.atom_count();
    Vector_3D_impl  zero (0,0,0);

    for (int i=0;  i<count;  ++i)
    {
        Atom & current_atom = atom_group_.atom(i);
        current_atom.kit().set_potential_force (zero);
    }
}

void Commonplace_interactions::
add_force ()
{
    int i, count;
    Atom_numerator numerator(atom_group_);

    try
    {
        count = atom_group_.atom_count();

        zero_forces ();      //fix remove?

        for (i=0;  i<count;  ++i)
        {
            Atom & current_atom = atom_group_.atom(i);
            add_atom_force (current_atom);
        }

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

            for (i=0;  i<count;  ++i)
            {
                Bond & current_bond = bond_group_.bond(i);
                add_bond_force (current_bond);
            }
        }
    }
    catch(...)
    {
        to_user().fatal_error ("Commonplace_interactions: Unknown exception.");
    }
}

void Commonplace_interactions::
atom_potential (Atom & atom)
{
    nonbonded_potential  (atom);
    angle_potential (atom);
}

void Commonplace_interactions::
bond_potential (Bond const & bond)
{
    bond_stretch_potential (bond);
    torsion_potential (bond);
}

void Commonplace_interactions::
add_atom_force (Atom & atom)
{
    add_nonbonded_force (atom);
    add_angle_force     (atom);
}

void Commonplace_interactions::
add_bond_force (Bond & bond)
{
    add_bond_stretch_force (bond);
    add_torsion_force      (bond);
}

void Commonplace_interactions::
nonbonded_potential (Atom & atom)
{
    double  lIepsilon = 1. / force_field().eps;
    double  A = 0., B = 0.;

    int     atom_number = atom.kit().number ();
    int     count       = atom_count();

    double  charge_332  = atom.charge() * EPS0; //331.8;
    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.kit().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 & current_atom  = atom_group_.atom (i);

        Force_field::VdWaals & tmp = vdw [current_atom.kit().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;

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

        P_el += current_atom.kit().transparency() *
            current_atom.charge() * charge_332 * lIr_eps;
    }

    process_14 (atom, false);
    restore_neighbours ();

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

int Commonplace_interactions::
atom_count ()
{
    return atom_group_.atom_count();
}

int Commonplace_interactions::
bond_count ()
{
    return bond_group_.bond_count();
}

}//MM

