#include "Pressure.h"

#include "Atom.h"
#include "Atom_kit.h"
#include "Interaction.h"
#include "Model_kit.h"

#include <math.h>

namespace MM
{

double Pressure::
value() const
{
    Interaction const& interaction = model_.kit().interaction();
    double virial = interaction.virial();

    int     count = model_.atom_count();
    double  sum = 0.;

    for (int i=0;  i<count;  ++i)
    {
        Atom const & atom = model_.atom (i);
        double mass = atom.mass();
        double velocity_2 = atom.kit().velocity().length_2();
        sum += mass * velocity_2;
    }

    // PV = N Kb T + <W>
    //      Tk
    double Tk = 0.5 * sum;

    Boundary_conditions const& boundary_conditions =
        model_.kit().boundary_conditions();

    if (!boundary_conditions.is_a_box())
        return 0;
    double volume = boundary_conditions.box().volume();

    double P = (virial + 2.*Tk) / (3.* volume);
    return P + correction();
}

double Pressure::
correction () const
{
    double Rcutoff = model_.kit().interaction().common().cutoff();
    return correction_ / (Rcutoff*Rcutoff*Rcutoff);
}

void   Pressure::
set_correction (double v)
{
    double Rcutoff = model_.kit().interaction().common().cutoff();
    correction_ = v * Rcutoff*Rcutoff*Rcutoff;
}

}//MM
