gnikit / md-sim

A mesh free molecular dynamics (MD) simulation for a fluid in an isothermal container

Molecular Dynamics Simulation

A program simulating Molecular Dynamics (MD) fluids, with the option to use custom pair potentials, with a range of initial lattice and boundary conditions.

Features

• Multiple pair potentials (easily extensible):
  • Bounded Inverse Power (BIP)
  • Gaussian Core Model (GCM)
  • Exponential Pair Potential (EXP)
  • Lennard-Jones Pair Potential (LJ)
• Multiple options for boundary conditions:
  • Periodic
  • Reflective
• Multiple initial lattice starting positions:
  • Simple Cubic (SC)
  • Face Centred Cubic (FCC)
  • Body Centred Cubic (BCC)
  • Random, based on normal distributions
• Multiple choices for iterative algorithms:
  • Explicit Verlet
  • Velocity Verlet
• Statistical quantities calculation:
  • Mean Square Displacement (MSD)
  • Velocity Autocorrelation Function (VAF)
  • Radial Distribution Function (RDF)
  • Structure Factor (SFx, SFy, SFz)
• 3D visualisation output through ParaView and CSV formatted logs
• Fluid compression for phase transitions
• Easy to use UI interface for setting up options
• Unit and regression testing

Getting Started

Clone

You can clone the project with all its submodules with:

git clone --recurse-submodules https://github.com/GNikit/md-sim.git

Build

Build the library and the md executable:

make -j

Testing

To run the test regression tests run:

make tests
make unit-tests

To add tests look into the tests folder.

Running a simulation

The preferred way for passing in options for a run is through the use of a .xml file. Such .xml files can be found in the tests directory or made from scratch using diamond (more on that later)

The following will run a simulation using an .xml input

./bin/md schemas/input_schema.xml

However, libmd can still by manually passing arguments as seen from the two examples below. Code without the use of the xml parser can be found in the examples directory.

Use libmd in your C++ project

A minimal working example of how to use the libmd.a is shown below, make sure to link and include the static library and the header file MD.h to your program.

#include "MD.h"

int main() {
  options_type options;

  options.steps = 5000;                     // number of iterations
  options.density = 0.5;                    // fluid density
  options.target_temperature = 0.5;         // fluid thermostat temperature
  options.potential_type = "LennardJones";  // pair potential type
  options.particles = {10, 10, 10};         // number of particles in xyz
  // More fine tuning options are available see data_structures.h

  MD run(options);
  run.Simulation();
}

Alternatively, one can use the non-options data structures for initialising MD

#include "MD.h"

int main() {
  size_t steps = 5000;  // number of iterations
  double rho = 0.5;     // fluid density
  double t = 0.5;       // fluid temperature

  MD run(steps, {10, 10, 10}, "SimpleCubic");
  run.Simulation("demo_", rho, t, NAN, NAN, "LennardJones");
}

Visualise using ParaView

To enable visualisation of the particles set the io/track_particles option to true in the schema, or directly set options.io_options.visualise = true before running a simulation.

Then load the xzy_data_...csv into ParaView, select Table of Points from the available filters, set x, y and z positions. For a nicer view, change the default representation Surface to Point Gaussian and increase the Gaussian Radius. A full set of instructions on how to load CSV files to Paraview can be found here.

Periodic Boundary conditions for a Lennard-Jones fluid	Reflective Boundary conditions for a Lennard-Jones fluid

Visualising log files

The log files produced from a run (if io options have been enabled) contain a range of statistical quantities indicative of how the fluid develops over time. The files are in CSV format, so almost any plotting tool can be used to visualise their values. A Python submodule in included under tools/md-tools that uses numpy and matplotlib, but the code is far from complete or robust.

Using diamond and libspud

This project relies on spud to easily read and set options from input files. Spud offers a visualisation tool called diamond to view the input options using a GUI. By default calling make on the top directory of md-sim compiles libspud and installs diamond and its helper functions in bin.

Use the following command to launch the GUI against a test_file.xml

./bin/diamond -s schemas/main_schema.rng test_file.xml

For more information about the diamond GUI see spud/doc/spud_manual.pdf

Updating the schema

If the schema files (rnc) are modified the rng files need to updated too. This can be done be calling:

make schema

md-sim's Issues

Add regression testing

Read output files and validate with runs that will compare their output with the saved data.
I was thinking of using a python interface rather than Catch2.

Particles spawn outside of box

Update readme

Update the info in README about:

• diamond and schemas
• remove the old mp4
• generate new gifs
• Update general info

Use more robust testing tools

The unittest module works but it is very annoying to add and debug test cases.
I was thinking of changing to pytest for regression testing
and Catch2 header only implementation for unit testing native c++.

Create new main.cpp to account for crystallisation

Add to the tinyxml2 interface in main.cpp that reads schema files and initialises a simulation
TODO list:

• Add crystallisation object in main.cpp
• Make tool generate crystallisation schemas

This will make it possible to launch all parts of the program from the command line

Make time step user defined

Allow for dt to be user defined, allowing/not allowing for temperature normalisation.

Temperature spike in step 0

A temperature spike is observed in the initial timestep, the calculated temperature is orders of magnitude greater than the target thermostat temperature.
The origin of such behaviour could indicate that the momentum of the system is not properly conserved. However, in the long run, the thermostat appears to take over and we reach the target temperature.

Infinite while loop for crystallization

It is possible to trigger an infinite while loop through the crystallisation method in phase_transition, by supplying a very large density increment

Change visualisation output to be compatible with ParaView

The idea is to use change to ASCII and binary CSV files written at every time step or at every collection of timesteps (user defined) so for the data files to be directly imported and visualised with ParaView.

This relates to issue #15

Add GUI support

Add a simple GUI, to handle the interactions between the vast number of options available

cutoff not resized during crystallisation

The cutoff radius is only assigned in the MD constructor hence when performing a compression run, the initial cutoff value is used, which in turn results into incorrect plots

Crystallisation segfault while outputting SF

SF is the wrong size when crystallising. Probably forgot to reset its size.

Change visualisation output to binary

When saving the visualisation coordinates in ASCII for python 3D animations can be time consuming and memory inefficient. Use a simple binary format for C++ or use some other binary fast storage method.

Add more stepping algorithms

RK4, etc.

Consider using PETSc

I am not entirely sure if this is worth the hassle, but I could use PETSc to store the position, velocity and force vectors in the form of (3D) Matrices.
This would allow

1. Easy parallelisation
2. Easy setup of a dynamic load balancer
3. Possible speed improvements

This efficient parallelisation could useful when tackling phase transitions.

Add melting step before crystallisation run

When attempting to run a crystallisation run, the temperatures are very low which most likely results into the initial lattice not melting and staying in the solid phase. Add the option to either

1. Melt the crystal with a high initial Temperature and then quench it.
2. Start with positions and velocities of a previous melted run.

Compression algorithm is faulty

I spotted the compression algorithm giving strange results, upon further inspection I saw that the box length was not scaled correctly.

Move to a more robust building system

It would be ideal if the project switched to cmake since currently the project uses handwritten Makefile which is far from ideal and prone to errors.