Vector 3 molecular dynamic with cell-list

Molecular Dynamic with Lennard-Jones potential

This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime. Particle feel each other by the potential.

\( V(x_p,x_q) = 4( (\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^6 ) \)

Constants

Here we define some useful constants

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

Calculate forces

In this function we calculate the forces between particles. It require the vector of particles Cell list and scaling factor for the Lennard-Jhones potential.

 
template<typename CellList> void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6, double r_cut2)
{
 

 
}
 

In the following we are going into the detail of this function

This function in called several time and require neighborhood of each particle. In order to speed-up the Cell-list construction we can use updateCellList function to reuse the memory of the previous cell-list. updateCellList can be faster than createCellList

See also

Cell-list

 
    vd.updateCellList(NN);
 

Get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q and calculate the force based on the Lennard-Jhones potential given by

\( F(x_p,x_q) = 24( \frac{2 \sigma^{12}}{r^{13}} - \frac{\sigma^6}{r^{7}}) \hat{r} \)

See also

Assign position

 
    // 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 force 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);
 
            if (rn > r_cut2)
            {++Np; continue;};
 
            // 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;
    }
 

Calculate energy

We also need a function to calculate energy, this function require the same parameter as calculate forces

 
template<typename CellList> double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList & NN, double sigma12, double sigma6, double r_cut2)
{
        double rc = r_cut2;
        double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
 

 
    // Calculated energy
    return E;
}
 

Reset the counter for the energy counter and update the cell list from the actual particle configuration

 
    double E = 0.0;
//  vd.updateCellList(NN);
 

In the following we are going into the detail of this function

First we get an iterator over the particles and get its position. For each particle p iterate in its neighborhood q and calculate the energy based on the Lennard-Jhones potential given by

\( V(x_p,x_q) = 4(\frac{1}{r^{12}} - \frac{1}{r^{6}}) \)

See also

Assign position

 
    // 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);
 
        // 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);
 
            if (rn > r_cut2)
            {++Np;continue;}
 
            // potential energy (using pow is slower)
            E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
 
            // 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;
    }
 

Initialization

After we defined the two main function calc forces and calc energy, we Initialize the library, we create a Box that define our domain, boundary conditions and ghost. Than we define important parameters of the simulation, time step integration, size of the box, and cut-off radius of the interaction. We also define 2 vectors x and y (they are like std::vector) used for statistic

See also

Initialization

 
    openfpm_init(&argc,&argv);
 
    double sigma = 0.1;
    double r_cut = 3.0*sigma;
 
    // we will use it do place particles on a 10x10x10 Grid like
    size_t sz[3] = {10,10,10};
 
    // domain
    Box<3,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
 
    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    // ghost, big enough to contain the interaction radius
    Ghost<3,double> ghost(r_cut);
 
    double dt = 0.0005;
    double sigma12 = pow(sigma,12);
    double sigma6 = pow(sigma,6);
 
    openfpm::vector<double> x;
    openfpm::vector<openfpm::vector<double>> y;
 

Than we define a distributed vector in 3D, containing 2 vectorial properties the first is the actual velocity of the particle the other is the force

See also

Vector instantiation

 
    vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost);
 

Particles on a grid like position

We define a grid iterator, to create particles on a grid like way. In the same cycle we also reset force and velocity

See also

Grid iterator

 
    // We create the grid iterator
    auto it = vd.getGridIterator(sz);
 
    while (it.isNext())
    {
        // Create a new particle
        vd.add();
 
        // key contain (i,j,k) index of the grid
        auto key = it.get();
 
        // The index of the grid can be accessed with key.get(0) == i, key.get(1) == j ...
        // We use getLastPos to set the position of the last particle added
        vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
        vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
        vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
 
        // We use getLastProp to set the property value of the last particle we added
        vd.template getLastProp<velocity>()[0] = 0.0;
        vd.template getLastProp<velocity>()[1] = 0.0;
        vd.template getLastProp<velocity>()[2] = 0.0;
 
        vd.template getLastProp<force>()[0] = 0.0;
        vd.template getLastProp<force>()[1] = 0.0;
        vd.template getLastProp<force>()[2] = 0.0;
 
        ++it;
    }
 
    vd.map();
    vd.ghost_get<>();
 

Molecular dynamic steps

Here we do 10000 MD steps using verlet integrator

The verlet integration stepping look like this

\[ \vec{v}(t_{n}+1/2) = \vec{v}_p(t_n) + \frac{1}{2} \delta t \vec{a}(t_n) \]

\[ \vec{x}(t_{n}+1) = \vec{x}_p(t_n) + \delta t \vec{v}(t_n+1/2) \]

calculate the forces from \( \vec{a} (t_{n}) \) finally

\[ \vec{v}(t_{n+1}) = \vec{v}_p(t_n+1/2) + \frac{1}{2} \delta t \vec{a}(t_n+1) \]

The cell-list structure is required to calculate forces. As demonstration purpose instead of using the standard Cell-list with (getCellList) we use the CELL_MEMBAL type. The impact in performance of using the CELL_MEMBAL instead of CELL_MEMFAST is less than 1% on the other hand CELL_MEMFAST use 16 Megabyte of memory

See also

Cell-list types

Inside this cycle we are using several features that has been explained before in particular

See also

Assign position

Mapping particles

Reduce (sum numbers across processors)

Ghost

Visualization, write VTK files

 
    // Get the Cell list structure
    auto NN = vd.getCellList<CELL_MEMBAL(3,double)>(r_cut);
 
    // The standard
    // auto NN = vd.getCellList(r_cut);
 
    // calculate forces
    calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
    unsigned long int f = 0;
 
    // MD time stepping
    for (size_t i = 0; i < 10000 ; i++)
    {
        // Get the iterator
        auto it3 = vd.getDomainIterator();
 
        // integrate velicity and space based on the calculated forces (Step1)
        while (it3.isNext())
        {
            auto p = it3.get();
 
            // here we calculate v(tn + 0.5)
            vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
            vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
            vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
 
            // here we calculate x(tn + 1)
            vd.getPos(p)[0] += vd.template getProp<velocity>(p)[0]*dt;
            vd.getPos(p)[1] += vd.template getProp<velocity>(p)[1]*dt;
            vd.getPos(p)[2] += vd.template getProp<velocity>(p)[2]*dt;
 
            ++it3;
        }
 
        // Because we moved the particles in space we have to map them and re-sync the ghost
        vd.map();
        vd.template ghost_get<>();
 
        // calculate forces or a(tn + 1) Step 2
        calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
 
        // Integrate the velocity Step 3
        auto it4 = vd.getDomainIterator();
 
        while (it4.isNext())
        {
            auto p = it4.get();
 
            // here we calculate v(tn + 1)
            vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
            vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
            vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
 
            ++it4;
        }
 
                // After every iteration collect some statistic about the configuration
                if (i % 100 == 0)
                {
                        // We write the particle position for visualization (Without ghost)
                        vd.deleteGhost();
                        vd.write_frame("particles_",f);
 
                        // we resync the ghost
                        vd.ghost_get<>();
 
                        // We calculate the energy
                        double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut*r_cut);
                        auto & vcl = create_vcluster();
                        vcl.sum(energy);
                        vcl.execute();
 
                        // we save the energy calculated at time step i c contain the time-step y contain the energy
                        x.add(i);
                        y.add({energy});
 
                        // We also print on terminal the value of the energy
                        // only one processor (master) write on terminal
                        if (vcl.getProcessUnitID() == 0)
                                std::cout << "Energy: " << energy << std::endl;
 
                        f++;
                }
    }
 

Plotting graphs

After we terminate the MD steps our vector x contains at which iteration we calculated the energy while y contains the energy value at that time-step. We can produce a graph X Y

Note

The graph produced is an svg graph that can be view with a browser. From the browser we can also easily save the graph into pure svg format

 
    // Google charts options, it store the options to draw the X Y graph
    GCoptions options;
 
    // Title of the graph
    options.title = std::string("Energy with time");
 
    // Y axis name
    options.yAxis = std::string("Energy");
 
    // X axis name
    options.xAxis = std::string("iteration");
 
    // width of the line
    options.lineWidth = 1.0;
 
    // Resolution in x
    options.width = 1280;
 
    // Resolution in y
    options.heigh = 720;
 
    // Add zoom capability
    options.more = GC_ZOOM;
 
    // Object that draw the X Y graph
    GoogleChart cg;
 
    // Add the graph
    // The graph that it produce is in svg format that can be opened on browser
    cg.AddLinesGraph(x,y,options);
 
    // Write into html format
    cg.write("gc_plot2_out.html");
 

Finalize

At the very end of the program we have always to de-initialize the library

 
    openfpm_finalize();
 

Full code

 
#include "Vector/vector_dist.hpp"
#include "Plot/GoogleChart.hpp"
#include "Plot/util.hpp"
#include "timer.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, double r_cut2)
{
 
 
 
    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 force 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);
 
            if (rn > r_cut2)
            {++Np; continue;};
 
            // 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 r_cut2)
{
        double rc = r_cut2;
        double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
 
 
 
    double E = 0.0;
//  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);
 
        // 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);
 
            if (rn > r_cut2)
            {++Np;continue;}
 
            // potential energy (using pow is slower)
            E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
 
            // 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;
    }
 
 
 
    // Calculated energy
    return E;
}
 
 
int main(int argc, char* argv[])
{
 
    openfpm_init(&argc,&argv);
 
    double sigma = 0.1;
    double r_cut = 3.0*sigma;
 
    // we will use it do place particles on a 10x10x10 Grid like
    size_t sz[3] = {10,10,10};
 
    // domain
    Box<3,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
 
    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    // ghost, big enough to contain the interaction radius
    Ghost<3,double> ghost(r_cut);
 
    double dt = 0.0005;
    double sigma12 = pow(sigma,12);
    double sigma6 = pow(sigma,6);
 
    openfpm::vector<double> x;
    openfpm::vector<openfpm::vector<double>> y;
 
 
 
    vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost);
 
 
 
    // We create the grid iterator
    auto it = vd.getGridIterator(sz);
 
    while (it.isNext())
    {
        // Create a new particle
        vd.add();
 
        // key contain (i,j,k) index of the grid
        auto key = it.get();
 
        // The index of the grid can be accessed with key.get(0) == i, key.get(1) == j ...
        // We use getLastPos to set the position of the last particle added
        vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
        vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
        vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
 
        // We use getLastProp to set the property value of the last particle we added
        vd.template getLastProp<velocity>()[0] = 0.0;
        vd.template getLastProp<velocity>()[1] = 0.0;
        vd.template getLastProp<velocity>()[2] = 0.0;
 
        vd.template getLastProp<force>()[0] = 0.0;
        vd.template getLastProp<force>()[1] = 0.0;
        vd.template getLastProp<force>()[2] = 0.0;
 
        ++it;
    }
 
    vd.map();
    vd.ghost_get<>();
 
 
    timer tsim;
    tsim.start();
 
 
    // Get the Cell list structure
    auto NN = vd.getCellList<CELL_MEMBAL(3,double)>(r_cut);
 
    // The standard
    // auto NN = vd.getCellList(r_cut);
 
    // calculate forces
    calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
    unsigned long int f = 0;
 
    // MD time stepping
    for (size_t i = 0; i < 10000 ; i++)
    {
        // Get the iterator
        auto it3 = vd.getDomainIterator();
 
        // integrate velicity and space based on the calculated forces (Step1)
        while (it3.isNext())
        {
            auto p = it3.get();
 
            // here we calculate v(tn + 0.5)
            vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
            vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
            vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
 
            // here we calculate x(tn + 1)
            vd.getPos(p)[0] += vd.template getProp<velocity>(p)[0]*dt;
            vd.getPos(p)[1] += vd.template getProp<velocity>(p)[1]*dt;
            vd.getPos(p)[2] += vd.template getProp<velocity>(p)[2]*dt;
 
            ++it3;
        }
 
        // Because we moved the particles in space we have to map them and re-sync the ghost
        vd.map();
        vd.template ghost_get<>();
 
        // calculate forces or a(tn + 1) Step 2
        calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
 
        // Integrate the velocity Step 3
        auto it4 = vd.getDomainIterator();
 
        while (it4.isNext())
        {
            auto p = it4.get();
 
            // here we calculate v(tn + 1)
            vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
            vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
            vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
 
            ++it4;
        }
 
                // After every iteration collect some statistic about the configuration
                if (i % 100 == 0)
                {
                        // We write the particle position for visualization (Without ghost)
                        vd.deleteGhost();
                        vd.write_frame("particles_",f);
 
                        // we resync the ghost
                        vd.ghost_get<>();
 
                        // We calculate the energy
                        double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut*r_cut);
                        auto & vcl = create_vcluster();
                        vcl.sum(energy);
                        vcl.execute();
 
                        // we save the energy calculated at time step i c contain the time-step y contain the energy
                        x.add(i);
                        y.add({energy});
 
                        // We also print on terminal the value of the energy
                        // only one processor (master) write on terminal
                        if (vcl.getProcessUnitID() == 0)
                                std::cout << "Energy: " << energy << std::endl;
 
                        f++;
                }
    }
 
 
    tsim.stop();
    std::cout << "Time: " << tsim.getwct() << std::endl;
 
 
    // Google charts options, it store the options to draw the X Y graph
    GCoptions options;
 
    // Title of the graph
    options.title = std::string("Energy with time");
 
    // Y axis name
    options.yAxis = std::string("Energy");
 
    // X axis name
    options.xAxis = std::string("iteration");
 
    // width of the line
    options.lineWidth = 1.0;
 
    // Resolution in x
    options.width = 1280;
 
    // Resolution in y
    options.heigh = 720;
 
    // Add zoom capability
    options.more = GC_ZOOM;
 
    // Object that draw the X Y graph
    GoogleChart cg;
 
    // Add the graph
    // The graph that it produce is in svg format that can be opened on browser
    cg.AddLinesGraph(x,y,options);
 
    // Write into html format
    cg.write("gc_plot2_out.html");
 
 
 
    openfpm_finalize();
 
 
}
 
 
 
 

Code with expression

Here we also show how we can simplify the example using expressions

See also

Vector 2 expression

#include "Vector/vector_dist.hpp"
#include "Plot/GoogleChart.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
#include "timer.hpp"
 
constexpr int velocity = 0;
constexpr int force = 1;
 
struct ln_potential
{
    double sigma12,sigma6,r_cut2,shift;
 
    ln_potential(double sigma12_, double sigma6_, double r_cut2_, double shift_) {sigma12 = sigma12_; sigma6 = sigma6_; r_cut2 = r_cut2_;shift = shift_;}
 
    Point<2,double> value(const Point<3,double> & xp, const Point<3,double> xq)
    {
        double rn = norm2(xp - xq);
        if (rn >= r_cut2)   return 0.0;
 
        Point<2,double> E({2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift,0.0});
 
        return E;
    }
};
 
struct ln_force
{
    double sigma12,sigma6,r_cut2;
 
    ln_force(double sigma12_, double sigma6_, double r_cut2_) {sigma12 = sigma12_; sigma6 = sigma6_;r_cut2 = r_cut2_;}
 
    Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)
    {
        Point<3,double> r = xp - xq;
        double rn = norm2(r);
 
        if (rn > r_cut2)    return 0.0;
 
        return 24.0*(2.0 * sigma12 / (rn*rn*rn*rn*rn*rn*rn) -  sigma6 / (rn*rn*rn*rn)) * r;
    }
};
 
int main(int argc, char* argv[])
{
    double dt = 0.0005, sigma = 0.1, r_cut = 3.0*sigma;
 
    double sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
    double rc2 = r_cut * r_cut;
    double shift = 2.0 * ( sigma12 / (rc2*rc2*rc2*rc2*rc2*rc2) - sigma6 / ( rc2*rc2*rc2) );
 
    openfpm::vector<double> x;
    openfpm::vector<openfpm::vector<double>> y;
 
 
    openfpm_init(&argc,&argv);
    Vcluster<> & vcl = create_vcluster();
 
    size_t sz[3] = {10,10,10};
    Box<3,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
    Ghost<3,double> ghost(r_cut);
 
    ln_force lf(sigma12,sigma6,r_cut*r_cut);
    ln_potential lp(sigma12,sigma6,r_cut*r_cut,shift);
 
    vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > vd(0,box,bc,ghost);
 
    auto v_force = getV<force>(vd);
    auto v_velocity = getV<velocity>(vd);
    auto v_pos = getV<PROP_POS>(vd);
 
    auto it = vd.getGridIterator(sz);
 
    while (it.isNext())
    {
        vd.add();
 
        auto key = it.get();
 
        vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
        vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
        vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
 
        ++it;
    }
 
    v_force = 0;
    v_velocity = 0;
 
    timer tsim;
    tsim.start();
 
    auto NN = vd.getCellList(r_cut);
 
    vd.updateCellList(NN);
    v_force = applyKernel_in_sim(vd,NN,lf);
    unsigned long int f = 0;
 
    // MD time stepping
    for (size_t i = 0; i < 10000 ; i++)
    {
        assign(v_velocity, v_velocity + 0.5*dt*v_force,
               v_pos, v_pos + v_velocity*dt);
 
        vd.map();
        vd.template ghost_get<>();
 
        // calculate forces or a(tn + 1) Step 2
        vd.updateCellList(NN);
        v_force = applyKernel_in_sim(vd,NN,lf);
 
        v_velocity = v_velocity + 0.5*dt*v_force;
 
        if (i % 100 == 0)
        {
            vd.deleteGhost();
            vd.write_frame("particles_",f);
            vd.ghost_get<>();
 
            vd.updateCellList(NN);
            Point<2,double> E = rsum(applyKernel_in_sim(vd,NN,lp) + (v_velocity * v_velocity)/2.0).get();
 
            vcl.sum(E.get(0));vcl.sum(E.get(1));
            vcl.execute();
 
            // we save the energy calculated at time step i c contain the time-step y contain the energy
            x.add(i);
            y.add({E.get(0),E.get(1),E.get(0) - E.get(1)});
 
            if (vcl.getProcessUnitID() == 0)
                std::cout << "Energy Total: " << E.get(0) << "   Kinetic: " <<  E.get(1) << "   Potential: " << E.get(0) - E.get(1) << std::endl;
 
            f++;
        }
    }
 
    tsim.stop();
    std::cout << "Time: " << tsim.getwct() << std::endl;
 
    GCoptions options;
    options.title = std::string("Energy with time");
    options.yAxis = std::string("Energy");
    options.xAxis = std::string("iteration");
    options.lineWidth = 1.0;
 
    GoogleChart cg;
    cg.AddLinesGraph(x,y,options);
    cg.write("gc_plot2_out.html");
 
    openfpm_finalize();
}