Introduction to molecular dynamics through OOP#
Molecular Dynamics allows to simulate complex systems where you can represent the interactions through forces and then solve the second Newton equation using some numerical method. It is applied on a large kind of problems, from cosmological simulations, to molecules and atoms, to biophysical systems, etc
|
|
In the previous section we learned how to solve this kind of initial value problems, using algorithms like rk4. But those algorithms require several evaluations of the force function per time step, and also become complex to implement when many particles are involved. So here we will better start changing the computational model to the object oriented approach and , while doing it, we will learn some concepts from the Molecular dynamics world. We will keep working on the particle bouncing on the floor problem. To do so, our simulation will have almost the same structure as before, but now we are talking about a particle:
Setup particle initial conditions
Iterate over time:
compute particle forces
perform a time step integration
print needed info
This looks similar to the previous approach but the mental model is simpler and also can be easily extended to more particles and/or complex forces. And to improve even more our code, we will separate the particle, from the integration, and from the forces. Each of these procedures will be on different files. Let’s start with the particle.
To learn more about classes, check https://www.learncpp.com/cpp-tutorial/classes-and-class-members/ .
Our goal now will be to be able to simulate a particle bouncing on the floor, with possibly more particles and walls. This could be easily applicable to more complex systems. This a basic physics engine, and we could later play to get something like https://www.youtube.com/watch?v=lS_qeBy3aQI&t=0s
%%html
<iframe width="560" height="315" src="https://www.youtube.com/embed/lS_qeBy3aQI" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
The particle class#
First of all, we will start thinking about the smallest modeling unit here, the
particle. To do so, we need to select the more fundamental characteristics of
our particle for the problem at hand, in this case, we need the particle to have
A mass, radius (type
double)A vector representing the position, the velocity and the forces. (type
std::vector<double>)
Many other characteristics could enter the simulation when the model needs them.
To declare this new data type, we use the following code by means of a struct
(we could also use a class, the only difference is that for structs all
members are public by default, while for a class all members are private by
default)
struct Particle {
std::vector<double> R{0, 0, 0}, V{0, 0, 0}, F{0, 0, 0};
double mass{0};
Particle(){
R = {0, 0, 0}; V = {0, 0, 0}; F = {0, 0, 0};
mass = 0.0;
}
};
here, R, V, and F, represent the position, velocity and force,
respectively. Each instance of the class Particle will have its own R, V
and so on. Let’s see this example in the python tutor or gdb online:
#include <vector>
struct Particle {
std::vector<double> R{0, 0, 0}, V{0, 0, 0}, F{0, 0, 0};
double mass{0};
Particle(){
R = {0, 0, 0}; V = {0, 0, 0}; F = {0, 0, 0};
mass = 0.0;
}
};
int main(int argc, char **argv)
{
Particle body1, body2;
body1.R[0] = 5.6;
body2.V[2] = -0.98;
return 0;
}
Adding methods#
This is only packing attributes, but you can also pack functions, like
#pragma once
#include <iostream>
#include <vector>
struct Particle {
std::vector<double> R{0,0,0}, V{0,0,0}, F{0,0,0};
double mass{0}, rad{0};
void print(void);
// overload the cout operator: friend declared to acces possible private data
// see: https://www.learncpp.com/cpp-tutorial/overloading-the-io-operators/
friend std::ostream& operator<< (std::ostream& out, const Particle & p) {
return out << p.mass << "\t" << p.rad << "\t"
<< p.R[0] << "\t" << p.R[1] << "\t" << p.R[2] << "\t"
<< p.V[0] << "\t" << p.V[1] << "\t" << p.V[2] << "\t"
<< p.F[0] << "\t" << p.F[1] << "\t" << p.F[2] << "\t";
}
};
For instance, there are also important member functions, like the constructor and the
destructor, that must be used when those tasks are important. You can also
overload operators to teach the language how to react, for instance, when you
try to sum up to particles (if that makes sense). In our case, the constructor
is already declared inside the class (Particle()), does not return anything
(not even void), and a destructor will be declared similarly but as
~Particle(). Using constructors/destructors is important when you need to
administer resources.
Our simulation will now be done using the particle struct. But our previous integration methods will need to be heavily adapted. It is much better for us to actually implement better integration algorithms right now, like the Leap-Frog algorithm.
Exercise#
Write the Particle struct into two files: Particle.h, with the needed
declarations, and Particle.cpp, with the implementations.
Implementing a better integration algorithm: leap-frog#
When studying ODE-IVP, we started with the Euler algorithm, which is pretty simple and computes the forces only once per time step but, as we have seen before, its order is low and it is not stable. In Molecular Dynamics there are several algorithms that fulfill several interesting properties (like being symplectic) that are commonly used, like leap-frog, verlet, velocity-verlet, optimized velocity verlet and so on. They try to integrate the evolution with a high order while reducing the number of force evaluations per time step. In particular, the leap-frog algorithm can be written as (https://en.wikipedia.org/wiki/Leapfrog_integration)
As you can see, the velocity is now in half-steps. This algorithm is very stable and symplectic. One important detail for the implementation: To perform the first step, you need to put your variables in the correct time. If we put t=0, the next step would be
Therefore, we will need first to move the velocity from time 0 to time \(0 - \delta t\). We could do that as follows
Exercise#
Create a new class called TimeIntegrator which receives on its constructor
the value of the time step, \(\delta t\). Implement there the Leap-frog algorithm
by using two member functions, startIntegration and timeStep, which receive
a array of particles (use templates for an arbitrary array of particles type,
particle_array_t) and perform the needed task. Implement everything inside the
class (only use Integrator.h)
#pragma once
class TimeIntegrator{
double dt{0.0};
public:
TimeIntegrator(double DT) { dt = DT; }
template <class particle_array_t>
void startIntegration(particle_array_t & parray) {
// this function moves the initial velocity to -dt/2
// V = V - A dt/2 = V - F dt /2m
for (auto & p : parray) {
for (int ii = 0; ii < 3; ++ii) {
p.V[ii] = p.V[ii] - p.F[ii]*dt/(2*p.mass);
}
}
}
template <class particle_array_t>
void timeStep(particle_array_t & parray) {
// TODO
}
};
Force computation#
Finally, we can start creating a function or object that computes the force for
a particle or for a set of particles. This normally called a collider, and it
will need some parameters (like the gravity constant, the floor stiffness and so
on) to perform its jobs. We can also use a simple function that receives the
parameters and the particle array and compute the force. Since the amount of
parameters is arbitrary, and their names, they can be stored in a dictionary
(std::map). They can also be split into particle level parameters (that can be
included on the particles themselves), or global parameters (like gravity). We
will here create a collider class and show how to implement the forces for our
simple case.
Exercise#
Create a new class called Collider which includes a member function to compute
external forces (for now) on an array of particles (arbitrary type, use
templates). To do so, complete the following code:
#pragma once
#include <map>
class Collider {
std::map<std::string, double> params; // parameters to compute the forces
public:
Collider(const std::map<std::string, double> & PARAMS) {
params = PARAMS;
}
template<class particle_array_t>
void computeForces(particle_array_t & parray) {
// reset forces
for (auto & p : parray) {
p.F = {0,0,0};
}
// individual forces
for (auto & p : parray) {
// gravity force
p.F[2] += p.mass*params["G"];
//TODO: force with floor
}
}
};
The main implementation#
Now we are ready to perform our simulation, using everything in our main function,
#include "Particle.h"
#include "Integrator.h"
#include "Collider.h"
#include <vector>
int main(int argc, char **argv) {
std::vector<Particle> bodies;
bodies.resize(1); // only one particle for now
// parameters
std::map<std::string, double> p;
p["T0"] = 0.0;
p["TF"] = 10.8767;
p["DT"] = 0.01;
p["G"] = -9.81;
// Force collider
Collider collider(p);
// Time initialization
TimeIntegrator integrator(p["DT"]);
// initial conditions and properties
bodies[0].R[2] = 0.987;
bodies[0].rad = 0.103;
bodies[0].mass = 0.337;
collider.computeForces(bodies); // force at t = 0
integrator.startIntegration(bodies); // start integration algorithm
std::cout << p["T0"] << "\t" << bodies[0] << "\n";
// Time iteration
const int niter = int((p["TF"] - p["T0"])/p["DT"]);
for(int ii = 1; ii < niter; ++ii) {
collider.computeForces(bodies);
integrator.timeStep(bodies);
double time = p["T0"] + ii*p["DT"];
std::cout << time << "\t" << bodies[0] << "\n";
}
return 0;
}
Exercises#
Add a damping with the ground. Where should you modify the code?
Plot the phase space with and without damping
Create an animation by printing the position data to a csv file and then using paraview.
Adding left and right walls#
Now put a wall on the left (-L/2) and on the right (+L/2). L is a new parameter. In the initial condition,s, make sure that you are NOT intersecting any wall at the beginning.
Add collisions with left and right walls. You will need to specify a length L. Put the left(right) wall at -L/2(L/2).
Adding more particles#
Now add another particle and implement a collision between particles. Later, Add many other particles and visualize their collisions. How does the cpu time scales with the number of particles?
If two particles are colliding, with interpenetration \(\delta\), then we can model the force that particle \(i\) exerts on particle \(j\) as,
where \(k\) is a constant proportional to the effective Young modulus, and \(\hat n_{ij}\) is the normal vector from particle \(i\) to \(j\), that is
with \(\vec R_{ij} = \vec R_j - \vec R_i\), and
For a more realistic force model, check https://www.cfdem.com/media/DEM/docu/gran_model_hertz.html and references therein.
Exercises#
Make the boundary a circle
Create a canon of particles to add particles to the system every nth step.
More about OOP modeling#
Modeling a complex number [7/10]#
Let’s model a complex number. The program should able to print the number and also to compute
[X] the real part.
[X] the imaginary part.
[X] the norm.
[X] the angle.
[X] the operator + .
[X] the operator - .
[ ] the operator / .
[ ] the operator * .
[X] the sin.
[ ] the cos.
Header file
#include <iostream>
#include <cmath>
struct Complex {
public:
double x{0.0}, y{0.0};
void print(void);
double real(void);
double imag(void);
double norm(void);
double arg(void);
Complex operator+(const Complex & rhs);
Complex operator-(const Complex & rhs);
Complex sin(void);
};
Implementation file
void Complex::print(void)
{
std::cout << "(" << x << " , " << y << ")" << std::endl;
}
double Complex::real(void)
{
return x;
}
double Complex::imag(void)
{
return y;
}
double Complex::norm(void)
{
return std::hypot(x, y);
}
double Complex::arg(void)
{
return std::atan2(y, x);
}
Complex Complex::operator+(const Complex & rhs)
{
Complex tmp;
tmp.x = x + rhs.x;
tmp.y = y + rhs.y;
return tmp;
}
Complex Complex::operator-(const Complex & rhs)
{
Complex tmp;
tmp.x = x - rhs.x;
tmp.y = y - rhs.y;
return tmp;
}
Complex Complex::sin(void)
{
Complex tmp;
tmp.x = std::sin(x)*std::cosh(y);
tmp.y = std::cos(x)*std::sinh(y);
return tmp;
}
Main file example
int main(void)
{
Complex z1, z2, z3;
z1.x = 4.32;
z1.y = 5.87;
z2.x = 3.0009;
z2.y = -2.65;
z1.print();
std::cout << "z1 real : " << z1.real() << std::endl;
std::cout << "z1 imag : " << z1.imag() << std::endl;
std::cout << "z1 norm : " << z1.norm() << std::endl;
std::cout << "z1 arg : " << z1.arg() << std::endl;
z1.print();
z2.print();
std::cout << "z3 = z1 + z2 : " << std::endl;
z3 = z1 + z2;
z3.print();
std::cout << "z3 = z1 - z2 : " << std::endl;
z3 = z1 - z2;
z3.print();
std::cout << "z3 = z1.sin() : " << std::endl;
z3 = z1.sin();
z3.print();
return 0;
}