DPD Exercise 1: DPD, hydrophobicity and parameterisation


DPD is particularly useful for looking at phase behaviour of multiple-component systems, because it evolves quickly in time and large-scale changes, such as phase separation, emerge quickly in the system. This is in contrast with molecular dynamics, where the relative hardness of particles greatly slows down the evolution of large-scale structures.

Groot and Warren’s [Groot1997] form of conservative interaction provides a quadratic potential and a quadratic equation of state - it is not particularly realistic but we can still parameterise it for single particle species by considering the isothermal compressibility of fluids.

We often want to use DPD to model how different species interact with each other. One approach we can use is to map the conservative force parameter between bead species to the Flory-Huggins solution theory of polymers. Assuming two particle species - A and B - interact among themselves in the same way, a parameter for the non-ideal part of the mixing Gibbs free energy between those species, \(\chi^{\text{AB}}\), can be related to the conservative interaction parameter between the two species by means of a proportionality:

(1)\[\chi^{\text{AB}} \propto \left( A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}} \right).\]

We can determine the proportionality constant from simulations of two segregating species of DPD beads. By keeping the overall bead density constant, varying the value of \(A_{ij}^{\text{AB}}\) and finding volume fractions of one species (e.g. \(\phi_{A}\) for species A) in bead-separated regions, we can use the following relationship to find \(\chi^{\text{AB}}\):

(2)\[\chi^{\text{AB}} = \frac{\ln\lbrack(1 - \phi_{A})/\phi_{A}\rbrack}{1 - 2\phi_{A}}.\]

Note that it may be difficult to obtain reliable results for mixed systems (i.e. when \(\chi^{\text{AB}} \lesssim 3\) or completely immiscible systems (i.e. when \(\chi^{\text{AB}}\) is rather large).

A relationship between \(\chi^{\text{AB}}\) and \(A_{ij}^{\text{AB}}\), will allow us to find conservative force parameters between pairs of species for which we know Flory-Huggins \(\chi\)-parameters. If we do not have experimental values of \(\chi\)-parameters available, it is possible to calculate them from energies of mixing determined from e.g. atomistic molecular dynamics simulations.

Please refer to Theory to DPD Exercise 1: DPD, hydrophobicity and parameterisation if you would like more details about the theoretical background to this Exercise.


We want to find the actual relationship between \(A_{ij}^{\text{AB}}\) and \(\chi^{\text{AB}}\) for DPD systems when \(\rho=3\). To do this, we need to carry out a series of calculations representing segregation between beads of species A and species B. From each of these calculations, we can find regions where one species is dominant over the other and work out the volume fraction of that species. We can then use that value in equation (2) to calculate a value of \(\chi^{\text{AB}}\) corresponding to the value of \(A_{ij}^{\text{AB}}\). By calculating \(\chi^{\text{AB}}\) for several values of \(A_{ij}^{\text{AB}}\), we can find a relationship between the two by fitting a line or curve using a regression technique.


Follow the instructions in Setting up DL_MESO to download DL_MESO and compile its DPD code (DL_MESO_DPD) and associated utilities. In particular for this exercise, you will need the main executable dpd.exe (either serial or parallel) as well as the utilities local.exe and either traject_vtf.exe or traject_xml.exe.

Copy the CONTROL, CONFIG and FIELD input files from the directory DEMO/DPD/FloryHuggins in your copy of DL_MESO to your working directory (e.g. dl_meso/WORK). These three input files represent a DPD simulation of a periodic cuboidal box with dimensions \(20 \times 8 \times 8\) and 3840 beads, half of which are of species A and the other half species B. (You might like to take a look at these files in a text editor and see how they are set out: Chapter 10 of the DL_MESO User Manual will give you more details about the formats.)

The CONTROL file will ask DL_MESO_DPD to run for 120,000 timesteps, of which the first 20,000 will be used to allow the system to settle into an equilibrated state before collecting trajectory data.

The CONFIG file initially positions the beads so that those of species A start in one half of the box and those of species B start in the other half. (We can run the simulation without this file, but DL_MESO_DPD will initially distribute the beads randomly throughout the box and it may therefore take longer for the system to reach our desired equilibrated state.)

The FIELD file specifies (among other things) values of \(A_{ij}\) between pairs of bead species: in this case, \(A_{ij}^{\text{AA}} = A_{ij}^{\text{BB}} = 25\) and \(A_{ij}^{\text{AB}} = 37\).

Follow the instructions in Setting up DL_MESO (the section ‘Running DL_MESO’) to run the simulation and then produce the traject.vtf file (or traject_*.xml files) using the traject_vtf.exe (or traject_xml.exe) utility.

To obtain the time-averaged plot of particle densities etc. by slicing up the simulation box, you can use the local.exe utility, which can be run as follows when slicing the box into 100 parts in the \(x\)-direction, 1 each in the \(y\)- and \(z\)-directions:

local.exe -nx 100 -ny 1 -nz 1 -av
./local.exe -nx 100 -ny 1 -nz 1 -av

which will produce a file called averages.vtk, that you can load into Paraview. (If you omit -av, you will also get the instantaneous density profiles for each recorded timeframe in the HISTORY file: these will be statistically very noisy!)

The later tasks can be carried out using all of the above, but you might find it easier and quicker to use the Python3 scripts, which are explained in more detail in Intro to Dissipative Particle Dynamics and can be downloaded here. The flory_huggins.py script will allow you to scan through ranges of conservative force parameters, obtain local concentration profiles and obtain values of \(\chi^{\text{AB}}\). The flory_huggins_plot.py script uses the results from the previous script and can plot the concentration profiles, recalculate values of \(\chi^{\text{AB}}\), plot the values of \(\chi^{\text{AB}}\) against \(A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}}\) and find the relationship.


  1. Run DL_MESO_DPD in your working directory with the supplied input files. Use the traject.exe utility to produce a traject.vtf file from the HISTORY file, then load the traject.vtf file into VMD. Take a look at the results: what happens to the different species?
  2. If the two bead species have separated, run the local.exe utility and split the box into several slices along the \(x\)-direction (e.g. 200) but only one each in the \(y\)- and \(z\)-directions.
  3. Open Paraview and load in the averages.vtk file (File->Open to load the file and then click Apply in Properties to create the visualisation).
    • Take a look at the density of one of the two species and produce a Surface plot. Can you see two distinct regions in the box?
    • Apply the Calculation filter to work out and visualise the total particle densities in each slice. Is the total density constant throughout the box? If there are any regions that are different to most of the box, what do you think is happening in those regions?
    • Apply the Calculation filter again to calculate the volume fraction of one species (i.e. find density_A/(density_A+density_B)) and then use the Plot Over Line filter along the \(x\)-axis to get a graph of this property.
    • Use the plot of volume fraction to work out a representative value of \(\phi_A\) for one of the bulk regions, then use equation (2) to calculate \(\chi^{\text{AB}}\).
  4. Run through a series of DL_MESO_DPD simulations, varying \(A_{ij}^{\text{AB}}\) from 33 to 43, obtain concentration profiles for each simulation, calculate \(\chi^{\text{AB}}\) values from those profiles and plot \(\chi^{\text{AB}}\) against \(\Delta A_{ij} = \left(A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}}\right)\). (Either do this by hand or use the flory_huggins.py Python script with its defaults and flory_huggins_plot.py to visualise both the concentration profiles and the mixing parameter data points.)
    • Do all of the simulations produce clearly separated regions for each particle species? Which ones do not manage this? (Why might that be?)
    • What kind of relationship exists between \(\Delta A_{ij} = \left(A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}}\right)\) and \(\chi^{\text{AB}}\)?
    • Is there a minimum value of \(\Delta A_{ij}\) that gives separation of the two species?
    • What relationships between \(\chi^{\text{AB}}\) and \(\Delta A_{ij}\) exist based on the simulations?
  5. Run some more simulations to calculate values of \(\chi^{\text{AB}}\), only this time either
    • Vary the value of \(A_{ij}^{\text{AA}}\) and the range of values for \(A_{ij}^{\text{AB}}\) to produce some more data points for the previous plots. (You may want to use a larger spacing between values of \(A_{ij}^{\text{AB}}\) to get through this more quickly.) Do these additional points still lie on more-or-less the same line as before?
    • Vary the system particle density from the default of 3.0 to a higher value (say, 5.0 or 6.0). Do you still get the same kind of relationship between \(\chi^{\text{AB}}\) and \(\Delta A_{ij}\) as before? How does the relationship change compared with the original particle density?