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
Apply constrains
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)
// particle.hpp
#pragma once
#include <iostream>
#include <valarray>
struct Particle {
std::valarray<double> R{0.0, 0.0, 0.0}, V{0.0, 0.0, 0.0}, F{0.0, 0.0, 0.0};
double mass{0}, rad{0};
void print(void);
};
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 "particle.h"
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
// particle.hpp
#pragma once
#include <iostream>
#include <valarray>
struct Particle {
std::valarray<double> R{0.0, 0.0, 0.0}, V{0.0, 0.0, 0.0}, F{0.0, 0.0, 0.0};
double mass{0}, rad{0};
void print(void);
};
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.
The actual implementartion will be in file particle.cpp
#include "particle.h"
// particle.cpp
void Particle::print(void) {
std::cout << mass << "\t" << rad << "\t"
<< R[0] << "\t" << R[1] << "\t" << R[2] << "\t"
<< V[0] << "\t" << V[1] << "\t" << V[2] << "\t"
<< F[0] << "\t" << F[1] << "\t" << F[2] << "\t";
}
// You can also overload the cout operator: friend declared to acces possible private data
// see: https://www.learncpp.com/cpp-tutorial/overloading-the-io-operators/
Notice the :: scope resolution operator.
Write the Particle struct into two files: Particle.h, with the needed
declarations, and Particle.cpp, with the implementations. Create a main file that can instantiate a particle and set some of its internal values, then print it.
Later just put everythin in the .hpp.
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
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, implement everything there)
#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) {
// this function moves the velocity by dt
// TODO
}
};
Constraints#
Normally, in this physical simulations we want to use physical forces for the system behaviour. But sometimes we don’t actually need a lot of precision, or to follow stric physical laws, so we can simplify a lot the interactions by applying “constraints”, i.e, forcing the particle to be in certain region, to be not mobile, etc.
In this case, we will play with keeping the particle inside a disc (sphere) of radius REXT, centering at the origin for simplicity. How to do that? at each step, we check the distance of the external border of our particle, and if it is larger than the external radius, we move the particle back to the inside of the region and invert its normal velocity. To do this, we will create a new class to amnage this kind of boundaries, called Boundary. It will have the needed parameters and will apply the needed constraints.
In our case, we want to keep all particles inside the sphere. To do so, for each particle, we can compute its external border with the spherical boundary, and detect if the particle is outside. Then, we can correct the particle position in the normal direction towards de center, and also invert the velocity, even adding a restitutin coefficient.
Deduce the vectorial corrections needed to apply the constraint described, to both the position and velocity, and implement them in the following class
#pragma once
#include <valarray>
#include <iostream>
class Boundary{
double RMAX_{0.0}, EN_{0.0}; // RMAX is the radius of the sphere, EN is the coefficient of restitution
std::valarray<double> C_{0.0, 0.0, 0.0}; // Center
public:
Boundary(double RMAX, double CX, double CY, double CZ, double EN) {
RMAX_ = RMAX;
C_[0] = CX;
C_[1] = CY;
C_[2] = CZ;
EN_ = EN;
}
template <class particle_array_t>
void apply(particle_array_t & parray) {
// applySphericalConstraint
for (auto & p : 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.
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,
// g++-14 -std=c++17 particle.cpp main_md.cpp
#include "particle.h"
#include "integrator.h"
#include "collider.h"
#include "boundary.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"] = 100.8767;
p["DT"] = 0.01;
p["G"] = 0.0; //-9.81;
// Force collider
Collider collider(p);
// Time initialization
TimeIntegrator integrator(p["DT"]);
// Boundary conditions
Boundary bc(2.345, 0.0, 0.0, 0.0, 1.0); // RMAX, CX, CY, CZ, EN
// initial conditions and properties
bodies[0].R[2] = 0.987; // z is upwards, x to the right
bodies[0].V[0] = 4.9876;//12.987; // z is upwards, x to the right
bodies[0].V[2] = 0.0; //4.9876; //3.987; // z is upwards, x to the right
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].print();
std::cout << "\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);
bc.apply(bodies);
double time = p["T0"] + ii*p["DT"];
std::cout << time << "\t";
bodies[0].print();
std::cout << "\n";
}
return 0;
}
Exercises#
Phase space#
Plot the phase space with and without damping
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?
Here we will implemement a physical force among particles. But, still, how would you do it using constraints?
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.
Add a particle cannon#
Create a cannon of particles to add particles to the system every nth step.
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).
A better visualization: using Paraview#
Paraview is “the world’s leading open source post-processing visualization engine. It integrates with your existing tools and workflows, allowing you to build visualizations to analyze data quickly. With its open, flexible, and intuitive user interface, you can analyze extremely large datasets interactively in 3D or programmatically using ParaView’s batch processing.” . Paraview can run on single laptops, or large scale clusters. You can use it to not only visualize but also analyse huge amounts of data. It has a large list of scientific data readers. Notice that there are more visualization tools, like visit, ovito, pymol, vmd, and so on. At the end it is important that you print your data in a format that can be read by these programs, like text, or better, something like vtk, hdf5, and so on.
Here, we will learn how to perform a very simple visualzation of our md data. We will need to :
Choose a format to print the data (text or binary? one file with all time steps or many files? which reader will we use? …)
Run the simulation and print the data
Use paraview to load the data, applying some filters to it.
Let’s do it.
Printing the data in a suitable format#
To keep this simple, we will add a funtion to print the system data. The filename will be the timestep, padded to , let’s say, 8 spaces. And the content, with the for each particle per line, would be
x, y, z, rad, mass, vx, vy, vz, fx, fy, fz
separated by , . Create the function and rerun the simulation.
NOTE: To avoid filling the hard disc and the folder with a lot of files, it is recommended to print the system info only every some steps, and in a separated folder. We will use a folder called DISPLAY. You can just create it in bash and assume its existence in the program, or use #include <filesystem>. and its utilities.
Loading the data into paraview#
Not that we have a lot of csv files inside the folder, let’s use paraview to open them. The path can change, but at the computer room you just run
/mnt/scratch/paraview/bin/paraview &
Then:
Load the data, with
openin the file many. Hit apply.Since this is a lot of csv files, we need to make them points with attributes. Go to filter and choose the filter
Table to Points. Select the right settings and press apply.Now, add a glyph, of sphere type, and as scale use the fild representing the radius. Now you have your system visualized and if you press play, you can see your system in action!
In paraview you can start applying filter, ike only printing particles with some specific speed, or adding some seed histogram in the right, or cut the filter in half and only plot the left part, and so on.
We can also perform an intermidiate step and use a python script to read our csv data and transform it into a vtk binary file, using pyevtk, and the load it in paraview. En sample script would be
import pandas as pd
import numpy as np
from pyevtk.hl import points_to_vtk
def csv_to_vtk(csv_filepath, vtk_filepath):
"""
Ingests a CSV file, reads particle data, and writes it to a VTK file
using pyevtk.
Args:
csv_filepath: Path to the input CSV file.
vtk_filepath: Path to the output VTK file.
"""
try:
# 1. Ingest CSV using pandas (handles potential variations in CSV format)
df = pd.read_csv(csv_filepath) # Let pandas infer data types
# 2. Extract data as NumPy arrays (required by pyevtk)
x = df['x'].to_numpy()
y = df['y'].to_numpy()
z = df['z'].to_numpy()
mass = df['mass'].to_numpy()
radius = df['radius'].to_numpy()
vx = df['vx'].to_numpy()
vy = df['vy'].to_numpy()
vz = df['vz'].to_numpy()
fx = df['fx'].to_numpy()
fy = df['fy'].to_numpy()
fz = df['fz'].to_numpy()
# 3. Write VTK using pyevtk
points_to_vtk(vtk_filepath, x, y, z,
data = {"mass": mass, "radius": radius,
"vx": vx, "vy": vy, "vz": vz,
"fx": fx, "fy": fy, "fz": fz})
print(f"Successfully wrote VTK file to: {vtk_filepath}")
except FileNotFoundError:
print(f"Error: CSV file not found at: {csv_filepath}")
except KeyError as e:
print(f"Error: Missing column in CSV: {e}")
except Exception as e: # Catch other potential errors
print(f"An error occurred: {e}")
# Example usage:
csv_file = "particle_data.csv" # Replace with your CSV file path
vtk_file = "particle_data.vtk" # Replace with desired VTK file path
csv_to_vtk(csv_file, vtk_file)