Molecular dynamic with symmetric interactions
In a previous example we show how to build a parallel molecular dynamic simulation. We also show how it was possible to achive better performance using Verlet list.
- See Also
- Molecular Dynamic with Lennard-Jones potential with verlet list
In this example we show how to improve even more our code using symmetric Verlet-list. If we look at the form of the Lennard-Jhones potential we will see that calculating the force that the particle p produce on q is equivalent to the force that q produce on p. This mean that when we use normal verlet-list we are redundantly doing double calculation. In order to avoid it we will use symmetric verlet-list. Symmetric verlet list have their feature in store the pair p,q only one time. If p store q as neighborhood than q does not store p as neighborhood. This Mean that when we calculate the contribution for of q to p we have to add such contribution also to q.
The example is exactly equivalent to the non-symmetric with few differences in calc_forces and calc_energies
Calculate forces
In this function we calculate the forces between particles. It require the vector of particles Cell list and sigma factor for the Lennard-Jhones potential. The function is exactly the same as the original
- See Also
- Calculate forces
with the following changes
- If we calculate the force for p-q we are also adding this force to q-p
vd.template getProp<force>(q)[0] -= f.get(0);
vd.template getProp<force>(q)[1] -= f.get(1);
vd.template getProp<force>(q)[2] -= f.get(2);
- At the end of the calculation we have to execute a ghost put
Explanation
The first point is given by the fact that if the pair is stored once, when we calculate the force, we have to add the contribution to both particles. The second instead is is given by the fact that q can be a ghost particles. In case q is a ghost particle we are adding the contribution to the ghost particle and not to the real one. To add the contribution to the real particle we have to use the function ghost_put. This function send back the information to the original processor, that will merge the information (in this case add_)
template<typename VerletList>
{
while (itg.isNext())
{
++itg;
}
while (it2.isNext())
{
auto p = it2.get();
Point<3,double> xp = vd.
getPos(p);
auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
while (Np.isNext())
{
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_cut * r_cut) {++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);
vd.template getProp<force>(q)[0] -= f.
get(0);
vd.template getProp<force>(q)[1] -= f.
get(1);
vd.template getProp<force>(q)[2] -= f.
get(2);
++Np;
}
++it2;
}
}
Calculate energy
For the energy we use symmetric verlet-list in the same way as we did for calc_forces. Because the symmetric verlet-list span each couple one time, we have to remove the division by two (in this case we use the original factor 4.0 of the Lennard-Jhones potential rather than 2.0).
template<typename VerletList>
{
double E = 0.0;
double rc = r_cut*r_cut;
double shift = 4.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
while (it2.isNext())
{
auto p = it2.get();
Point<3,double> xp = vd.
getPos(p);
auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
double Ep = E;
while (Np.isNext())
{
if (q == p.getKey()) {++Np; continue;};
Point<3,double> xq = vd.
getPos(q);
double rn = norm2(xp - xq);
if (rn >= r_cut*r_cut) {++Np;continue;}
E += 4.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;
}
Simulation
The simulation is equal to the simulation explained in the example molecular dynamic
- See Also
- Molecular Dynamic with Lennard-Jones potential with verlet list
The difference is that we create a symmetric Verlet-list instead of a normal one
The rest of the code remain unchanged
double dt = 0.00025;
double sigma = 0.1;
double r_cut = 3.0*sigma;
double r_gskin = 1.3*r_cut;
double sigma12 = pow(sigma,12);
double sigma6 = pow(sigma,6);
openfpm_init(&argc,&argv);
size_t sz[3] = {10,10,10};
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
size_t k = 0;
size_t start = vd.
accum();
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;
k++;
++it;
}
calc_forces(vd,NN,sigma12,sigma6,r_cut);
unsigned long int f = 0;
int cnt = 0;
double max_disp = 0.0;
for (size_t i = 0; i < 10000 ; i++)
{
double max_displ = 0.0;
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];
Point<3,double> disp({vd.template getProp<velocity>(p)[0]*dt,vd.template getProp<velocity>(p)[1]*dt,vd.template getProp<velocity>(p)[2]*dt});
vd.
getPos(p)[0] += disp.get(0);
vd.
getPos(p)[1] += disp.get(1);
vd.
getPos(p)[2] += disp.get(2);
if (disp.norm() > max_displ)
max_displ = disp.norm();
++it3;
}
if (max_disp < max_displ)
max_disp = max_displ;
if (cnt % 10 == 0)
{
vd.template ghost_get<>();
}
else
{
vd.template ghost_get<>(SKIP_LABELLING);
}
cnt++;
calc_forces(vd,NN,sigma12,sigma6,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);
auto & vcl = create_vcluster();
x.add(i);
y.add({energy});
std::cout << "Energy: " << energy << " " << max_disp << " " << std::endl;
max_disp = 0.0;
f++;
}
}
std::cout <<
"Time: " << tsim.
getwct() << std::endl;
Finalize
At the very end of the program we have always to de-initialize the library
Full code
#include "Vector/vector_dist.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "data_type/aggregate.hpp"
#include "Plot/GoogleChart.hpp"
#include "Plot/util.hpp"
#include "timer.hpp"
constexpr int velocity = 0;
constexpr int force = 1;
template<typename VerletList>
{
while (itg.isNext())
{
++itg;
}
while (it2.isNext())
{
auto p = it2.get();
Point<3,double> xp = vd.
getPos(p);
auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
while (Np.isNext())
{
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_cut * r_cut) {++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);
vd.template getProp<force>(q)[0] -= f.
get(0);
vd.template getProp<force>(q)[1] -= f.
get(1);
vd.template getProp<force>(q)[2] -= f.
get(2);
++Np;
}
++it2;
}
}
template<typename VerletList>
{
double E = 0.0;
double rc = r_cut*r_cut;
double shift = 4.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
while (it2.isNext())
{
auto p = it2.get();
Point<3,double> xp = vd.
getPos(p);
auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
double Ep = E;
while (Np.isNext())
{
if (q == p.getKey()) {++Np; continue;};
Point<3,double> xq = vd.
getPos(q);
double rn = norm2(xp - xq);
if (rn >= r_cut*r_cut) {++Np;continue;}
E += 4.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[])
{
double dt = 0.00025;
double sigma = 0.1;
double r_cut = 3.0*sigma;
double r_gskin = 1.3*r_cut;
double sigma12 = pow(sigma,12);
double sigma6 = pow(sigma,6);
openfpm_init(&argc,&argv);
size_t sz[3] = {10,10,10};
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
size_t k = 0;
size_t start = vd.
accum();
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;
k++;
++it;
}
calc_forces(vd,NN,sigma12,sigma6,r_cut);
unsigned long int f = 0;
int cnt = 0;
double max_disp = 0.0;
for (size_t i = 0; i < 10000 ; i++)
{
double max_displ = 0.0;
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];
Point<3,double> disp({vd.template getProp<velocity>(p)[0]*dt,vd.template getProp<velocity>(p)[1]*dt,vd.template getProp<velocity>(p)[2]*dt});
vd.
getPos(p)[0] += disp.get(0);
vd.
getPos(p)[1] += disp.get(1);
vd.
getPos(p)[2] += disp.get(2);
if (disp.norm() > max_displ)
max_displ = disp.norm();
++it3;
}
if (max_disp < max_displ)
max_disp = max_displ;
if (cnt % 10 == 0)
{
vd.template ghost_get<>();
}
else
{
vd.template ghost_get<>(SKIP_LABELLING);
}
cnt++;
calc_forces(vd,NN,sigma12,sigma6,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);
auto & vcl = create_vcluster();
x.add(i);
y.add({energy});
std::cout << "Energy: " << energy << " " << max_disp << " " << std::endl;
max_disp = 0.0;
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();
}