DPD Tutorial Exercise 1: Simple DPD systems¶

Ex2. Parameterising DPD¶

Most DPD simulations involve multiple components and how they interact with each other, form structures at equilibrium, react to flow fluids etc. If we want to use DPD to model a given system, we need to understand the effect of conservative interactions between different particle (bead) species and how we can put this to use. In order to do this, we now need to start using DL_MESO_DPD.

We have supplied a ZIP file with the source code for both of DL_MESO’s codes. Unpack this ZIP file using the command:

[ ]:

!unzip -o -q dl_meso_2.7.zip

This produces a set of directories, including the DL_MESO_DPD source code in dl_meso/DPD and a working directory of its own in dl_meso/WORK. To compile the serial (one processor core) version of DL_MESO_DPD, use the following commands to navigate into the dl_meso/WORK directory and invoke the makefile:

[ ]:

%%bash
cd dl_meso/WORK
make -f Makefile-serial clean
make -f Makefile-serial

This provides an executable file for DL_MESO_DPD (dpd.exe) that we can use to run a simulation, provided we have the required input files.

For our first DL_MESO_DPD simulation, we have a CONTROL and a FIELD file available in a working directory (DPD1Ex2). To launch the calculation, either use the following commands (ideally in a terminal window):

[ ]:

%%bash
cd DPD1Ex2
../dl_meso/WORK/dpd.exe

or the following script to launch DL_MESO_DPD in the specified folder and keep track of its progress.

[ ]:

import launchdlmeso as dlm
dlm.run_DPD('DPD1Ex2', 'dl_meso/WORK/dpd.exe', 1, True, '')

DL_MESO_DPD will take a couple of minutes to run and produce four output files:

CORREL - a statistics file with values of propeties at set intervals during the calculation after equilibration
HISTORY - the simulation trajectory with simulation snapshots (configurations), written in binary format for speed and intended to be processed using utilities or scripts
export - simulation restart file with the current simulation state, used to resume the simulation at a later date or to start a new simulation from the given state
REVIVE - simulation restart file with statistical accumulators and random number generator states, used to resume the simulation at a later date

N.B. DL_MESO_DPD normally prints information about the calculation (e.g. timing data, instantaneous system properties) while it runs to an OUTPUT file. To get the script with the progress bar to work properly we have added the l_scr option to the CONTROL file to redirect that output to the screen or standard output, which is then directed (or ‘piped’) to a file called OUTPUT by the script.

To plot the data in the CORREL file, either use plotting software or a spreadsheet program, or the following plotting script.

[ ]:

from plotcorrel import *
plotCorrel('DPD1Ex2/CORREL', 'time', 'pe-total', 'Time', 'Potential energy per particle', 'Potential energy per particle vs. time')
plotCorrel('DPD1Ex2/CORREL', 'time', 'temperature', 'Time', 'Temperature', 'Temperature vs. time')
plotCorrel('DPD1Ex2/CORREL', 'time', 'pressure', 'Time', 'Pressure', 'Pressure vs. time')

To visualise the simulation, we can use the HISTORY file to generate a file (traject.vtf) that can be read into VMD. We can do this by using the following script:

[ ]:

%run history_dlm_to_vtf.py --in DPD1Ex2/HISTORY --out DPD1Ex2/traject.vtf

Once you have managed to run DL_MESO_DPD, you might want to consider the following:

If you look in the FIELD file, you will see that \(A_{ij}^{AB}\), the conservative force parameter between beads of species A and B, is signficantly higher than those for each species (\(A_{ij}^{AA}\) and \(A_{ij}^{BB}\)). Given DL_MESO_DPD initially mixes up beads of both species in the box (look at the first trajectory frame), what happens to them during the simulation?
Try gradually reducing \(A_{ij}^{AB}\) and re-run the calculation. What happens when this value gets close to \(A_{ij}^{AA}\)? When it dips below \(A_{ij}^{AA}\)?

Going on … we want to see how the value of \(A_{ij}^{AB}\) affects separation of two particle species. To do this, we can carry out a series of simulations that use different values of the conservative force parameter, measure the volume fractions of one species (e.g. \(\phi_A\)) in separated regions (away from interfaces where the species mix), calculate the Flory-Huggins parameter:

\[\chi^{AB} = \frac{\ln [(1-\phi_A)/\phi_A]}{1 - 2\phi_A}\]

and plot a graph of \(A_{ij}^{AB}\) vs. \(\chi^{AB}\).

Rather than doing this manually, we have come up with a Python script that will help automate this process. Before we go ahead and launch it, you might want to re-compile DL_MESO_DPD to use OpenMP multithreading to speed up each DPD calculation: you can do this by using the following commands with the required Makefile.

[ ]:

%%bash
cd dl_meso/WORK
make -f Makefile-OMP clean
make -f Makefile-OMP

Use the following commands to launch the workflow script. It will launch a series of DL_MESO_DPD calculations in an output folder using different values of \(A_{ij}^{AB}\) and work out corresponding values of \(\chi^{AB}\) from time-averaged profiles of concentrations (volume fractions). The script will run 11 calculations, each of which will take a few minutes, so the whole process might take a little while to complete. (Hence why we suggested re-compiling DL_MESO_DPD with OpenMP!)

Note that the number of OpenMP threads you want to use for each calculation is set using the environment variable OMP_NUM_THREADS using the command just before the script is launched. We do not recommend using all available threads (16 in this working environment) as there are overheads involved in setting them up and using them, but 4 threads per calculation should work well enough.

(Optionally: If you want to see how we have set up the calculations, take a look at the CONFIG file generated by the script in VMD. You will need to open it as a ‘DL_POLY CONFIG’ file: note that both DL_POLY and DL_MESO_DPD use the same file format.)

[ ]:

import os
os.environ["OMP_NUM_THREADS"] = "4"
%run flory_huggins.py --dlmeso dl_meso/WORK/dpd.exe

The flory_huggins.py script will produce a text file with each simulation’s concentration profile (\(\phi_A (x)\) vs. \(x\)), the values of \(A_{ij}^{AA}\) and \(A_{ij}^{AB}\) used to generate it and the value of \(\chi^{AB}\) with an error estimate. To plot each of the concentration profiles, invoke the following plotting script (changing the last number to select each profile).

[ ]:

from plotfloryhuggins import *
plotConcentration('out/floryhuggins-rho-3.000.dat', 5)

If you just cannot wait for these results, we have provided a sample output file based on using the flory_huggins.py script: floryhuggins-rho-3.000.dat. You can modify the location of the output file in the above and below calls to the plotting script (i.e. remove out/) to use this file instead of the one you end up generating.

Quick question based on these results:

Can you see clearly separated regions for all of the simulations? If not, which simulations did not manage this and why might this be?

Another script is available to plot the obtained values of \(\chi^{AB}\) against \(\Delta A = A_{ij}^{AB} - A_{ij}^{AA}\) and find a best-fit line that goes through the origin (i.e. \(\chi^{AB} = 0\) when \(A_{ij}^{AB} = A_{ij}^{AA}\)), adding the equations for this line as a title.

[ ]:

from plotfloryhuggins import *
plotChi('out/floryhuggins-rho-3.000.dat')

Questions and tasks based on this result:

Groot and Warren in their 1997 Journal of Chemical Physics suggested that \(\chi^{AB} \propto \Delta A\) … do we get this?
We can add more concentration profiles and data points to the above plot by re-running our Python script, only this time try changing the value of \(A_{ij}^{AA}\) and the ranges of \(A_{ij}^{AB}\) values, e.g. for \(A_{ij}^{AA} = 50.0\), \(A_{ij}^{AB}\) between 58 and 68, and spacing out the \(A_{ij}^{AB}\) values by 2.0 rather than 1.0:

[ ]:

import os
os.environ["OMP_NUM_THREADS"] = "4"
%run flory_huggins.py --dlmeso dl_meso/WORK/dpd.exe --Aii 50.0 --Aijmin 58.0 --Aijmax 68.0 --dA 2.0

Does the value of \(A_{ij}^{AA}\) have an effect on the relationship between \(\chi^{AB}\) and \(\Delta A\)?
(Optional) We have so far assumed an overall particle density \(\rho = 3\), but we can try a different density using another command-line option for the flory_huggins.py script, --rho, e.g. the below command for \(\rho = 5\). (Note the different range of \(A_{ij}^{AB}\) values and that this will write the results to a different file, so you will need to change the filenames in the above plotConcentration and plotChi function calls.)

[ ]:

import os
os.environ["OMP_NUM_THREADS"] = "4"
%run flory_huggins.py --dlmeso dl_meso/WORK/dpd.exe --rho 5.0 --Aijmin 29.0 --Aijmax 36.0

What happens to the relationship between \(\chi^{AB}\) and \(\Delta A\) when \(\rho\) is changed?
Finally … how can we make use of these relationships to set up a DPD simulation?