#include "Harmonic_potential.h"

#include "Atom.h"
#include "Atom_kit.h"
#include "Point_3D.h"

namespace MM
{

Harmonic_potential::
Harmonic_potential (Atom & atom1, Atom & atom2)
:   
    atom1_(atom1), 
    atom2_(atom2),
    length_(atom1.position(), atom2.position()),
    force_constant_(1), 
    equilibrium_distance_(0)
{
}

Harmonic_potential::
Harmonic_potential (Atom & atom1,                 Atom & atom2,
                    double equilibrium_distance,  double force_constant)
:
    atom1_(atom1),
    atom2_(atom2),
    length_(atom1.position(), atom2.position()),
    force_constant_(force_constant),
    equilibrium_distance_(equilibrium_distance)
{
}

double Harmonic_potential::
potential ()
{
    double x = length_.value() - equilibrium_distance_;
    last_potential_ = force_constant_ * x * x;
    return last_potential_;
}

void Harmonic_potential::
add_force ()
{
    Vector_3D_impl direction_12 (0,0,0);
    atom2_.position().difference (atom1_.position(), &direction_12);
    double distance = direction_12.length();
    double f = -2. * force_constant_ * (distance - equilibrium_distance_);
    direction_12 *= f / distance;
    atom2_.kit().add_potential_force (direction_12);
    atom1_.kit().add_potential_force (direction_12.negate());
}

}//MM

