Using DL_MESO_DPD

Summary

DL_MESO_DPD uses a number of input files to read in the details of a DPD simulation and produces output files containing the results. The Java GUI and utilities supplied with DL_MESO can be used to create/modify DPD simulation inputs and work with outputs to enable users to visualise and analyse simulations. Three Python3 scripts are also supplied to make it easier to carry out a simulation workflow for one of the DPD practical exercises and examine the contents of an output file produced by DL_MESO_DPD.

Input and output files

The required input files for DL_MESO_DPD are:

  • CONTROL: gives the basic simulation controls, specifying properties such as timestep size, temperature, number of timesteps, ensemble type, boundary conditions etc.
  • FIELD: provides definitions and numbers of beads, as well as interactions between them, the contents and connectivity of molecules and any external fields.

These files are similar in form to the equivalent files in DL_POLY, with some small differences to make running DPD simulations easier. These can be created or modified either by hand (using your favourite text editor) or by using the Java GUI.

DL_MESO_DPD will devise an initial configuration (particle positions, velocities and forces) based on the contents of the FIELD file. If you wish to use a different initial configuration, you can also supply a CONFIG file, which is formatted very similarly to those DL_POLY. This file can be created from a previous simulation using some of the utilities supplied with DL_MESO_DPD.

The output files from DL_MESO_DPD can include:

  • OUTPUT: gives a summary of the DPD simulation, including input information, timing, statistical properties (e.g. energies, pressure, temperature) at given timesteps and averages over the simulation.
  • CORREL: gives tables of statistical properties at user-specified time intervals - can be read and plotted using e.g. Excel, gnuplot.
  • HISTORY: trajectory data (i.e. simulation snapshots at user-specified time intervals) - written in binary, requires utilities to read and use.
  • export: state of simulation recorded at time intervals, continuously overwritten (last one provided) - written in binary, can be used to restart simulation or used with utilities.
  • REVIVE: statistical accumulators and random number generator states for simulation, continuously overwritten - written in binary, can be used to restart simulation.

The output file formats for DL_MESO_DPD (HISTORY, export and CORREL) are deliberately agnostic, in the sense that they each do not generate data to be read by a single program. The utilities supplied with DL_MESO (see below) can convert this data into specific formats, but if you are able to write your own programs or scripts, you can get the data into the format of your choice.

Further details of these file formats can be found in Chapter 10 of the DL_MESO User Manual.

DL_MESO_DPD utilities

Most of the DPD utilities will read the HISTORY or export files and produce more useable outputs. The most useful ones are:

  • traject_vtf.exe: reads the HISTORY file’s trajectories and writes a trajectory file (traject.vtf) that can be read and visualised by VMD;
  • traject_xml.exe: reads the HISTORY file’s trajectories and writes a series of trajectory files (traject_*.xml) that can be read and visualised by OVITO;
  • local.exe: calculates instantaneous and time-averaged localised properties (densities, velocities, stress tensors) from HISTORY file trajectories by summing properties in user-specified slices of the system (in a similar fashion to histograms), can be visualised and analysed using Paraview;
  • isosurfaces.exe: calculates isosurfaces of a specified particle species (that can be visualised in Paraview) and an order parameter to determine mesophases;
  • rdf.exe: calculates radial distribution functions (RDFs) from HISTORY file trajectories;
  • export_config.exe: converts the simulation snapshot in export files into a CONFIG file (to allow new simulations starting from that configuration);
  • molecule.exe: determines inputs required to include molecules in DPD simulations, including sample configurations, bond connectivity etc. as required in the FIELD file.

Analysis and visualisation tools

For the exercises supplied, we recommend the following analysis and visualisation programs:

  1. VMD - this is a molecular visualisation program that can display, animate and analyse large particle-based systems. [1]
  2. OVITO - this is another molecular visualisation program that can display, animate and analyse large particle-based systems.
  3. Paraview - this is a visualisation program that can display, animate and analyse large grid-based datasets.
  4. Gnuplot - this is a flexible command-line based tool that can generate graphs from tabulated data.

All three are available for all operating systems and can either be downloaded directly from their websites or installed using package managers.

Python scripts

To help with the DPD exercises, I have provided three Python3 scripts for you to use. These can help you automate simulation workflows, visualise and analyse the results of DL_MESO_DPD simulations: they are not absolutely essential and you could still do the exercises without them, but they will speed things up!

These scripts use the following Python3 modules that you may need to install:

  • docopt
  • tqdm
  • numpy
  • PyQt5
  • matplotlib

These can be installed using the command pip3 install [module] [2]. Two of the Python scripts will also require the graphical user interface toolkit Qt 5, which is available for any operating system.

If you would like to make use of these scripts, we suggest putting a copy of each one into your working directory with the various executables for DL_MESO_DPD.

flory_huggins.py

flory_huggins.py is a simulation workflow script that will carry out a series of DL_MESO_DPD simulations to obtain a relationship between DPD conservative force parameters and the Flory-Huggins \(\chi\)-parameters (see DPD Exercise 1: DPD, hydrophobicity and parameterisation for more details).

It creates the required input files - CONTROL for simulation parameters, CONFIG for the initial particle configuration and FIELD for interaction parameters - to carry out these simulations for various values of the conservative force parameter between two different bead species, launches the simulations one-by-one, obtains a time-averaged local concentration profile for one bead species per simulation and adds these profiles to a text file for later analysis. It also calculates a value of \(\chi\) for each simulation and estimates its error.

To use the script, you can either invoke it with your Python3 interpreter:

python3 flory_huggins.py

or, if the permissions on the script allow it to be executed, you can use (in Unix/Linux or macOS):

./flory_huggins.py

You can add a number of optional command-line options at the end of the invocation:

  • --rho <rho> sets the overall particle density of the simulations to <rho> (equal to 3 by default)
  • --Aii <Aii> sets the conservative force parameter between beads of the same species to <Aii> (equal to 25.0 by default)
  • --Aijmin <Aijmin> sets the minimum conservative force parameter between beads of the two different species to <Aijmin> (equal to 33.0 by default)
  • --Aijmax <Aijmax> sets the maximum conservative force parameter between beads of the two different species to <Aijmax> (equal to 43.0 by default)
  • --dA <dA> sets the interval of unlike species conservative force parameter between simulations to <dA> (equal to 1.0 by default)
  • --dx <dx> sets the spacing of histogram bins for concentration profiles to <dx> (equal to 0.1 by default)
  • --L <L> sets the length of the box along the x-direction to <L> (equal to 20.0 by default)
  • --W <W> sets the width of the box in the y- and z-directions to <W> (equal to 8.0 by default)
  • --print-to-screen will redirect the summary of each simulation to the screen instead of writing to an OUTPUT file

The results will be placed in a file called floryhuggins-rho-<rho>.dat. The results of each simulation will consist of a header line containing the conservative force parameters for like-like and unlike bead interactions, an estimate of \(\chi\) with a simple measure of its standard deviation and the number of lines in the concentration profile that will follow, which is tabulated with \(x\) and \(\phi (x)\). Note that while different system densities will be put into separate files, results for different values of like-like conservative parameter will be collected together: subsequent runs of the script will append further data to this file. (If you want to start over from scratch, you will need to delete the file.)

The script currently assumes you wish to use DL_MESO_DPD on a single core for each simulation: if you have compiled the code with MPI and want to change this, you can modify the command in line 164. (See Setting up DL_MESO for more information.)

flory_huggins_plot.py

flory_huggins_plot.py will visualise the results from flory_huggins.py using the generated text file. You can look at the local concentration profiles generated from the DL_MESO_DPD simulations, recalculate \(\chi\) values for each set of conservative force parameters, plot these values and find a relationship between \(\chi^{AB}\) and \(A_{ij}^{AB} - A_{ij}^{AA}\).

To use the script, you can either invoke it with your Python3 interpreter:

python3 flory_huggins_plot.py

or, if the permissions on the script allow it to be executed, you can use:

./flory_huggins_plot.py

You can optionally add the command-line option --rho <rho>, which chooses data for the system-wide particle density <rho> (which is equal to 3.0 by default).

image0
Concentration plot

The Concentration tab allows you to plot the local concentration profile for each simulation: a pulldown box at the bottom left of the window will allow you to choose the simulation. The legend box includes the value of \(\chi^{AB}\) obtained for this simulation, along with an estimate of the standard deviation of this value. (This has not been determined in a very systematic way, but does give some measure of how accurate the value of \(\chi^{AB}\) happens to be.)

If you want to modify this value, you can choose to recalculate it by choosing the extents used to sample the concentrations: two regions are supplied (one for high concentration, one for low concentration) as indicated by vertical dashed lines. The positions of these lines can be changed using the text boxes at the bottom of the window, and the recalculate button can be used to obtain a new value (and standard deviation) for \(\chi^{AB}\). The reset button can be used to go back to the original values of \(\chi^{AB}\) as provided with the concentration data.

image1
Mixing parameter plot

The Mixing parameter tab gives a plot of \(\chi^{AB}\) as a function of \(\Delta A_{ij} = A_{ij}^{AB} - A_{ij}^{AA}\). If you have used multiple values of the like-like conservative force parameter \(A_{ij}^{AA}\), data for different values of this parameter will be plotted in different colours (as indicated by the legend box in the top left). The standard deviations for each data point are indicated using error bars.

Since the relationship between \(\chi^{AB}\) and \(\Delta A_{ij}\) is assumed to be entirely proportional, a dashed line indicating the best fit of the available data to a linear plot with zero at the origin is included. The equations for this best-fit line are given at the top of the plot, which include error bars indicating the extremes of the gradient.

Since the data below and above certain values of \(\Delta A_{ij}\) may not be reliable, you can change the minimum and maximum values of this property used to calculate the best-fit line in the text boxes at the bottom of the window. Changing these may modify both the line’s gradient and the equations shown at the top of the graph.

correl.py

correl.py is a more general-purpose script that will plot and analyse the results from DL_MESO_DPD, using its CORREL and (if available) RDFDAT files with statistical information and radial distribution functions respectively. (This is based on a similar script written for DL_POLY STATIS files called statis.py.)

To use the script, you can either invoke it with your Python3 interpreter:

python3 correl.py

or, if the permissions on the script allow it to be executed, you can use:

./correl.py

You can add up to two optional command-line options at the end of the invocation:

  • --cor <cor> sets the name of the CORREL file to <cor> (default: CORREL)
  • --rdf <rdf> sets the name of the RDFDAT file to <rdf> (default: RDFDAT)
image2
CORREL file plot

When the script is launched, the data in the CORREL file is read and is made available for plotting and analysis in the first tab. There are two pulldown boxes at the bottom left of the window: Data set selects which statistical property you wish to use - the list of which depends on the contents of the CORREL file - while Analysis selects what you want to do with the selected data. You can plot one of the following for a given statistical property:

  • A simple Timeseries
  • A Histogram of the values obtained during the simulation
  • A time series with a Running average added (also providing a measure of mean and standard error)
  • An Autocorrelation with respect to time
  • A Fourier transform
image3
RDFDAT file plot

If an RDFDAT file is available (by running the rdf.exe utility after a DL_MESO_DPD simulation), the second tab for this file will be populated with a plotting window. The Data set pulldown box allows you to select which pair of particle species you wish to plot, while the Analysis pulldown box allows you to choose to plot the raw Radial Distribution Function or its Fourier Transform (related to the structure factor). The plot for the latter can also depend upon the number of bins used to calculate the Fourier transform, which can be modified using the text box at the bottom of the window.

Footnotes

[1]If you are using the most recent version of macOS (Catalina 10.15 or newer), you will currently not be able to run VMD - at least, not without some kind of workaround - as it is a 32-bit application that relies on OpenGL.
[2]Depending on your Python3 installation, you may need root (administrator) permissions to use this command. In Unix/Linux and macOS, you can often do this by either inserting sudo before the command itself, which will prompt you for your user password, or by typing su to login as root (which requires the root password).