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
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
while (it2.isNext())
{
auto p = it2.get();
Point<3,double> xp = vd.
getPos(p);
vd.template getProp<force>(p)[0] = 0.0;
vd.template getProp<force>(p)[1] = 0.0;
vd.template getProp<force>(p)[2] = 0.0;
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.
getPos(p)));
while (Np.isNext())
{
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
Point<3,double> xq = vd.
getPos(q);
Point<3,double> r = xp - xq;
double rn = norm2(r);
if (rn > r_cut2)
{++Np; continue;};
Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
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);
++Np;
}
++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) );
Reset the counter for the energy counter and update the cell list from the actual particle configuration
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
while (it2.isNext())
{
auto p = it2.get();
Point<3,double> xp = vd.
getPos(p);
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.
getPos(p)));
while (Np.isNext())
{
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
Point<3,double> xq = vd.
getPos(q);
double rn = norm2(xp - xq);
if (rn > r_cut2)
{++Np;continue;}
E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
++Np;
}
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;
++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;
size_t sz[3] = {10,10,10};
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
double dt = 0.0005;
double sigma12 = pow(sigma,12);
double sigma6 = pow(sigma,6);
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
while (it.isNext())
{
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;
}
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
auto NN = vd.
getCellList<CELL_MEMBAL(3,
double)>(r_cut);
calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
unsigned long int f = 0;
for (size_t i = 0; i < 10000 ; i++)
{
while (it3.isNext())
{
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];
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;
}
vd.template ghost_get<>();
calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
while (it4.isNext())
{
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;
}
if (i % 100 == 0)
{
vd.
write(
"particles_",f);
double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut*r_cut);
auto & vcl = create_vcluster();
x.add(i);
y.add({energy});
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
options.
title = std::string(
"Energy with time");
options.
yAxis = std::string(
"Energy");
options.
xAxis = std::string(
"iteration");
cg.
write(
"gc_plot2_out.html");
Finalize
At the very end of the program we have always to de-initialize the library
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)
{
while (it2.isNext())
{
auto p = it2.get();
Point<3,double> xp = vd.
getPos(p);
vd.template getProp<force>(p)[0] = 0.0;
vd.template getProp<force>(p)[1] = 0.0;
vd.template getProp<force>(p)[2] = 0.0;
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.
getPos(p)));
while (Np.isNext())
{
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
Point<3,double> xq = vd.
getPos(q);
Point<3,double> r = xp - xq;
double rn = norm2(r);
if (rn > r_cut2)
{++Np; continue;};
Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
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);
++Np;
}
++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;
while (it2.isNext())
{
auto p = it2.get();
Point<3,double> xp = vd.
getPos(p);
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.
getPos(p)));
while (Np.isNext())
{
auto q = Np.get();
if (q == p.getKey()) {++Np; continue;};
Point<3,double> xq = vd.
getPos(q);
double rn = norm2(xp - xq);
if (rn > r_cut2)
{++Np;continue;}
E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
++Np;
}
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;
++it2;
}
return E;
}
int main(int argc, char* argv[])
{
openfpm_init(&argc,&argv);
double sigma = 0.1;
double r_cut = 3.0*sigma;
size_t sz[3] = {10,10,10};
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
double dt = 0.0005;
double sigma12 = pow(sigma,12);
double sigma6 = pow(sigma,6);
while (it.isNext())
{
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;
}
auto NN = vd.
getCellList<CELL_MEMBAL(3,
double)>(r_cut);
calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
unsigned long int f = 0;
for (size_t i = 0; i < 10000 ; i++)
{
while (it3.isNext())
{
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];
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;
}
vd.template ghost_get<>();
calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
while (it4.isNext())
{
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;
}
if (i % 100 == 0)
{
vd.
write(
"particles_",f);
double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut*r_cut);
auto & vcl = create_vcluster();
x.add(i);
y.add({energy});
std::cout << "Energy: " << energy << std::endl;
f++;
}
}
std::cout <<
"Time: " << tsim.
getwct() << std::endl;
options.
title = std::string(
"Energy with time");
options.
yAxis = std::string(
"Energy");
options.
xAxis = std::string(
"iteration");
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;
{
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;
}
};
{
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);
size_t sz[3] = {10,10,10};
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
ln_force lf(sigma12,sigma6,r_cut*r_cut);
auto v_force = getV<force>(vd);
auto v_velocity = getV<velocity>(vd);
auto v_pos = getV<PROP_POS>(vd);
while (it.isNext())
{
++it;
}
v_force = 0;
v_velocity = 0;
v_force = applyKernel_in_sim(vd,NN,lf);
unsigned long int f = 0;
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.template ghost_get<>();
v_force = applyKernel_in_sim(vd,NN,lf);
v_velocity = v_velocity + 0.5*dt*v_force;
if (i % 100 == 0)
{
vd.
write(
"particles_",f);
Point<2,double> E = rsum(applyKernel_in_sim(vd,NN,lp) + (v_velocity * v_velocity)/2.0,vd).get();
x.add(i);
std::cout <<
"Energy Total: " << E.
get(0) <<
" Kinetic: " << E.
get(1) <<
" Potential: " << E.
get(0) - E.
get(1) << std::endl;
f++;
}
}
std::cout <<
"Time: " << tsim.
getwct() << std::endl;
options.
title = std::string(
"Energy with time");
options.
yAxis = std::string(
"Energy");
options.
xAxis = std::string(
"iteration");
cg.
write(
"gc_plot2_out.html");
openfpm_finalize();
}