Sarkas: A Fast Pure-Python Molecular Dynamics Code

Sarkas is a pure-python molecular dynamics (MD) code we are currently developing for wide (open source) use. To our knowledge, Sarkas is the first production-scale (using Numba) pure-Python molecular dynamics code for simulating plasmas. Our current areas of application are strongly coupled plasmas; strongly coupled plasmas include dusty plasmas, ultracold neutral plasmas and warm dense matter.  Sarkas offers the ease of use of Python while being highly performant with execution speeds comparable to or exceeding that of compiled languages (e.g., C). Sarkas’s high-performance originates with extensive use of Numpy arrays and Numba’s just-in-time compilation.
Here, we show how MD can be used to study kinetic instabilities in plasmas. In this movie, the vertical axis is particle velocity and the horizontal axis is space (that is, it is phase space). The plasma is mostly composed of particles in equilibrium (red points) with a bump on the tail (blue particles at large, positive velocity). This situation mimics a plasma with an injected pulse of fast particles in it. What is interesting about this MD simulation is that it resolves particles down to their atomic scale, while successfully generating longer wavelength instabilities, thereby spanning extreme length scales. This shows promise for using Sarkas to study a wider range of kinetic and hydrodynamic phenomena building from the atomistic level.
Interparticle interactions currently included are long-range Coulomb and intermediate- to short-range screened Coulomb (Yukawa interactions). The computational cell is 3D periodic cubic.  For short range potentials, O(N) scaling is achieved with a neighbor list (linked cell); conversely, for intermediate- to long-range forces  the Particle-Particle Particle-Mesh (PPPM) algorithm is employed.  Details of the theoretical and numerical aspects can be found in the paper linked below.
As an example, MD simulations of strongly coupled plasmas enable access to the plasma’s static and dynamic properties. Static structural information is obtained through the radial distribution function g(r) using the particle positions. Dynamical properties such as transport processes  are also available by computing the transport coefficients from the particle trajectories. As a example, in the figure below the results of dynamic structure factor S(k,ω) for a Yukawa plasma, and the resulting wave dispersion properties.
skw_sarkas_web
To illustrate performance of Sarkas we compared the execution times with a equivalent version of the code written entirely in C. The figure below compares C and Python timings of two computationally intensive parts of the code: LCL for the Particle-Particle part and charge assignment to a grid for the Particle-Mesh part. The times(Intel Xeon E5-2670v2, 2.6GHz) are for one step of the force calculation for 104 – 106 particles. The performance of Sarkas is similar to that of the C version, even with C’s O3 compiler optimization. This is a remarkable result considering that Sarkas is a pure Python (with Numba) code.
pp_pm_sarkas_optC_ss_2
You can learn more about our recent work here:
Screen Shot 2017-03-25 at 5.38.17 PM
In its current form Sarkas uses Ewald generalized to intermediate range interactions. Many extensions to the code are expected. Stay tuned!

This work is lead by Gautham Dharuman.

Publications

A generalized Ewald decomposition for screened Coulomb interactions
Journal Chemical Physics. 146, 024112 (2017)
Gautham Dharuman, Liam G. Stanton, James N. Glosli and Michael S. Murillo