energy_force.hpp

/*
 * energy_force.hpp
 *
 *  Created on: Aug 6, 2016
 *      Author: i-bird
 */

#ifndef EXAMPLE_VECTOR_4_REORDER_ENERGY_FORCE_HPP_
#define EXAMPLE_VECTOR_4_REORDER_ENERGY_FORCE_HPP_

constexpr int velocity = 0;
constexpr int force = 1;


template<typename CellList> void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6)
{
    // update the Cell-list
    vd.updateCellList(NN);

    // Get an iterator over particles
    auto it2 = vd.getDomainIterator();

    // For each particle p ...
    while (it2.isNext())
    {
        // ... get the particle p
        auto p = it2.get();

        // Get the position xp of the particle
        Point<3,double> xp = vd.getPos(p);

        // Reset the forice counter
        vd.template getProp<force>(p)[0] = 0.0;
        vd.template getProp<force>(p)[1] = 0.0;
        vd.template getProp<force>(p)[2] = 0.0;

        // Get an iterator over the neighborhood particles of p
        auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));

        // For each neighborhood particle ...
        while (Np.isNext())
        {
            // ... q
            auto q = Np.get();

            // if (p == q) skip this particle
            if (q == p.getKey())    {++Np; continue;};

            // Get the position of p
            Point<3,double> xq = vd.getPos(q);

            // Get the distance between p and q
            Point<3,double> r = xp - xq;

            // take the norm of this vector
            double rn = norm2(r);

            // Calculate the force, using pow is slower
            Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) -  sigma6 / (rn*rn*rn*rn)) * r;

            // we sum the force produced by q on p
            vd.template getProp<force>(p)[0] += f.get(0);
            vd.template getProp<force>(p)[1] += f.get(1);
            vd.template getProp<force>(p)[2] += f.get(2);

            // Next neighborhood
            ++Np;
        }

        // Next particle
        ++it2;
    }
}


template<typename CellList> double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6)
{
    double E = 0.0;

    // Update the cell-list
    vd.updateCellList(NN);

    // Get the iterator
    auto it2 = vd.getDomainIterator();

    // For each particle ...
    while (it2.isNext())
    {
        // ... p
        auto p = it2.get();

        // Get the position of the particle p
        Point<3,double> xp = vd.getPos(p);

        // Reset the force
        vd.template getProp<force>(p)[0] = 0.0;
        vd.template getProp<force>(p)[1] = 0.0;
        vd.template getProp<force>(p)[2] = 0.0;

        // Get an iterator over the neighborhood of the particle p
        auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));

        // For each neighborhood of the particle p
        while (Np.isNext())
        {
            // Neighborhood particle q
            auto q = Np.get();

            // if p == q skip this particle
            if (q == p.getKey())    {++Np; continue;};

            // Get position of the particle q
            Point<3,double> xq = vd.getPos(q);

            // take the normalized direction
            double rn = norm2(xp - xq);

            // potential energy (using pow is slower)
            E += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) );

            // Next neighborhood
            ++Np;
        }

        // Kinetic energy of the particle given by its actual speed
        E +=   (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
                vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
                vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;

        // Next Particle
        ++it2;
    }

    return E;
}


#endif /* EXAMPLE_VECTOR_4_REORDER_ENERGY_FORCE_HPP_ */