OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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)
{
Class for FAST cell list implementation.
Definition CellList.hpp:357
Distributed vector.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
}

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;
}
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28

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);
This class represent an N-dimensional box.
Definition Box.hpp:61
Implementation of 1-D std::vector like structure.

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

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
// 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");
Small class to produce graph with Google chart in HTML.
void write(std::string file)
It write the graphs on file in html format using Google charts.
void AddLinesGraph(openfpm::vector< X > &x, openfpm::vector< Y > &y, const GCoptions &opt)
Add a simple lines graph.
Google chart options.
size_t width
width of the graph in pixels
size_t heigh
height of the graph in pixels
std::string xAxis
X axis name.
size_t lineWidth
Width of the line.
std::string more
more
std::string title
Title of the chart.
std::string yAxis
Y axis name.

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);
// 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
// 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();
}
Class for cpu time benchmarking.
Definition timer.hpp:28
void stop()
Stop the timer.
Definition timer.hpp:119
void start()
Start the timer.
Definition timer.hpp:90
double getwct()
Return the elapsed real time.
Definition timer.hpp:130

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;
{
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_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);
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;
cg.AddLinesGraph(x,y,options);
cg.write("gc_plot2_out.html");
openfpm_finalize();
}
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition Point.hpp:172
void execute()
Execute all the requests.
void sum(T &num)
Sum the numbers across all processors and get the result.
size_t getProcessUnitID()
Get the process unit id.
Implementation of VCluster class.
Definition VCluster.hpp:59