From 4d4dd9d1209a02c6fee2c1df29a6a76d80f889c6 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Sat, 14 Apr 2018 22:41:46 +0200
Subject: [PATCH] Adding missing file
---
example/Vector/3_molecular_dynamic/config.cfg | 2 +-
.../3_molecular_dynamic/main_expr_paper.cpp | 107 ++++++++++++++++++
2 files changed, 108 insertions(+), 1 deletion(-)
create mode 100644 example/Vector/3_molecular_dynamic/main_expr_paper.cpp
diff --git a/example/Vector/3_molecular_dynamic/config.cfg b/example/Vector/3_molecular_dynamic/config.cfg
index affb9522..5853db22 100644
--- a/example/Vector/3_molecular_dynamic/config.cfg
+++ b/example/Vector/3_molecular_dynamic/config.cfg
@@ -1,2 +1,2 @@
[pack]
-files = main.cpp main_expr.cpp main_vl.cpp Makefile
+files = main.cpp main_expr.cpp main_expr_paper.cpp main_vl.cpp Makefile
diff --git a/example/Vector/3_molecular_dynamic/main_expr_paper.cpp b/example/Vector/3_molecular_dynamic/main_expr_paper.cpp
new file mode 100644
index 00000000..cfe7b19e
--- /dev/null
+++ b/example/Vector/3_molecular_dynamic/main_expr_paper.cpp
@@ -0,0 +1,107 @@
+#include "Vector/vector_dist.hpp"
+#include "Plot/GoogleChart.hpp"
+#include "Operators/Vector/vector_dist_operators.hpp"
+#include "timer.hpp"
+
+
+void Init_grid(size_t (& sz)[3],
+ vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > & particles)
+{
+ auto it = particles.getGridIterator(sz);
+
+ while (it.isNext())
+ {
+ particles.add();
+
+ auto key = it.get();
+
+ particles.getLastPos()[0] = key.get(0) * it.getSpacing(0);
+ particles.getLastPos()[1] = key.get(1) * it.getSpacing(1);
+ particles.getLastPos()[2] = key.get(2) * it.getSpacing(2);
+
+ ++it;
+ }
+}
+
+double sigma12;
+double sigma6;
+double r_cut2;
+
+constexpr int velocity_prop = 0;
+constexpr int force_prop = 1;
+
+///// Define lennard_jones interaction to be used in applyKernel_in_sim
+DEFINE_INTERACTION_3D(ln_force)
+ 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;
+END_INTERACTION
+
+int main(int argc, char* argv[])
+{
+ ///// Initialize the library /////
+ openfpm_init(&argc,&argv);
+
+ ///// Initialize constants and initial configuration /////
+ double dt = 0.0005, sigma = 0.1, r_cut = 3.0*sigma;
+ sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
+ double rc2 = r_cut * r_cut;
+
+ ///// Define initialization grid, simulation box, periodicity
+ ///// and ghost
+ size_t sz[3] = {10,10,10};
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+ Ghost<3,float> ghost(r_cut);
+
+ ///// Lennard_jones potential object used in applyKernel_in
+ ln_force lennard_jones;
+
+ ////// Define particle and initialize on a grid
+ ////// the value 0 mean we start with 0 particles, Init_grid
+ ////// create them
+ vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > particles(0,box,bc,ghost);
+ Init_grid(sz,particles);
+
+ // Take aliases for the particle force, velocity and position
+ // we use them to write math expressions
+ auto force = getV<force_prop>(particles);
+ auto velocity = getV<velocity_prop>(particles);
+ auto position = getV<PROP_POS>(particles);
+
+ // set the velocity to zero (a math expression)
+ velocity = 0;
+
+ ///// Get the cell-list and calculate the force using lennard_jones
+ ///// potential
+ auto NN = particles.getCellListSym(r_cut);
+ force = applyKernel_in_sim(particles,NN,lennard_jones);
+
+ // Molecula Dynamic time stepping
+ for (size_t i = 0; i < 10000 ; i++)
+ {
+ // 1-step Verlet velocity
+ // v(t + 1/2*dt) = v(t) + 1/2*force(t)*dt
+ // x(t + dt) = x(t) + v(t + 1/2*dt)
+ velocity = velocity + 0.5*dt*force;
+ position = position + velocity*dt;
+
+ // redistribute particles and fill ghost
+ particles.map();
+ particles.ghost_get<>();
+
+ // Calculate the force at t + dt
+ particles.updateCellListSym(NN);
+ force = applyKernel_in_sim(particles,NN,lennard_jones);
+
+ // 2-step Verlet velocity
+ // v(t+dt) = v(t + 1/2*dt) + 1/2*force(t+dt)*dt
+ velocity = velocity + 0.5*dt*force;
+ }
+
+ openfpm_finalize();
+}
+
--
GitLab