DPD Exercise 1: DPD, hydrophobicity and parameterisation

Introduction

As a mesoscale modelling method, DPD is particularly useful for looking at phase and structural behaviour of multiple-component systems due to the method’s ability to step quickly in time. Large-scale changes such as phase separation and self-assembly can emerge more rapidly than they normally can with molecular dynamics, both in terms of available timescales and computational time.

A large part of this comes down to the choice of interaction potential. Groot and Warren’s [Groot1997] form of conservative interaction provides a quadratic potential and a quadratic equation of state. For single component systems, while this interaction limits us to a single fluid phase, it does incorporate some non-ideality and we can parameterise it by considering the fluid’s isothermal compressibility.

To extend DPD and Groot-Warren interactions to multiple components, we need to find some way to obtain conservative force parameters between pairs of different species. To that end, one approach available to us is to map these parameters to the Flory-Huggins solution theory of polymers.

If we assume particle species A and B interact among themselves in the same way (i.e. \(A_{ij}^{\text{AA}} = A_{ij}^{\text{BB}}\)), we can define a parameter for the non-ideal part of the mixing Gibbs free energy between A and B, \(\chi^{\text{AB}}\). This parameter can be shown to be proportional to an increase in conservative force parameter from the value for like-like interactions, i.e.

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

To determine the proportionality constant, we can carry out DPD simulations of two separating species of beads with different values of \(A_{ij}^{\text{AB}}\) and measure the resulting values of \(\chi^{\text{AB}}\). If the two species separate into regions rich in either species A or species B, we can find an average volume fraction of one species (e.g. \(\phi_{\text{A}}\) for species A) away from the interfaces between the regions. The volume fraction can then be used to calculate \(\chi^{\text{AB}}\):

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

Once we have the actual relationship between \(\chi^{\text{AB}}\) and \(A_{ij}^{\text{AB}}\), we can then calculate conservative force parameters between pairs of species for which we know lory-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.

More theoretical details on parameterising DPD simulations using Flory-Huggins solution theory are available as additional background for this Exercise.

Aim

We want to find the actual relationship between \(A_{ij}^{\text{AB}}\) and \(\chi^{\text{AB}}\) for DPD systems with a commonly-used particle density, \(\rho=3\).

To do this, we can follow the procedure given above to calculate \(\chi^{\text{AB}}\) from simulations of separating beads with different values of \(A_{ij}^{\text{AB}}\), holding both the overall bead density and the like-like conservative force parameters (\(A_{ij}^{\text{AA}}\) and \(A_{ij}^{\text{BB}}\)) constant. Using several data points, we can plot \(\chi^{\text{AB}}\) as a function of \(\left( A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}} \right)\) and get a best-fit line to find the proportionality constant for the overall bead density.

Instructions

If you have not yet done so, 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. Optionally, you can also make use of the flory_huggins.py and flory_huggins_plot.py Python scripts. You will also need to have Paraview and either VMD or OVITO installed.

Copy the CONTROL, CONFIG and FIELD input files from the directory dl_meso/DEMO/DPD/FloryHuggins in your copy of DL_MESO to your working directory (dl_meso/WORK). You might like to take a look at these files in a text editor and see how they are set out: Using DL_MESO_DPD has some more details on how CONTROL and FIELD are formatted.

In summary, these files represent a DPD simulation of a periodic cuboidal box with \(20 \times 8 \times 8\) and 3840 beads, half of which are of species A and the other half species B. We are using the CONFIG file to override DL_MESO_DPD’s default randomised setup and initially place the beads of species A in one half of the box and the beads for species B in the other half. We will be running DL_MESO_DPD for 120,000 timesteps, using the first 20,000 to let the system equilibrate before collecting trajectory data in a HISTORY file every 100 timesteps. The FIELD file specifies the interactions between beads: in this case, \(A_{ij}^{\text{AA}} = A_{ij}^{\text{BB}} = 25\) and \(A_{ij}^{\text{AB}} = 37\). (We are also setting \(r_c = 1\) for the interactions and \(\gamma = 1\) for all species pairs.)

To run the simulation, either type in:

./dpd.exe

if you have compiled the serial version of DL_MESO_DPD (with or without OpenMP), or if you want to run the parallel version on X processor cores, type:

mpirun -np X ./dpd.exe

Once the simulation has finished, type in either one of the following two commands:

./traject_vtf.exe
./traject_xml.exe

Both utilities will convert the contents of the HISTORY file into a format that can be read using visualisation software: traject_vtf.exe will generate a traject.vtf file to open in VMD, while traject_xml.exe will generate a series of files (traject_*.xml) that can be read by OVITO. (Which one to use will depend on what you have installed!)

A useful analysis we can carry out for this simulation is to slice the simulation box into a number of sections along the longest dimension (in this case, \(x\)), count the numbers of particles for each species and divide by the volume of each slice to get local particle densities. We can get these by running the local.exe utility using the following command:

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

This utility requires information about the number of slices in each dimension, which can be provided by the -nx, -ny and -nz command line options. (The above will split the box into 200 slices along the \(x\)-direction, but not split it in the \(y\)- or \(z\)-directions.) The -av command line option also asks local.exe to only provide a time-averaged profile in a file called averages.vtk: omitting this will also produce instantaneous profiles, although these will be statistically very noisy!

Once you have obtained the averages.vtk file, you can open it using Paraview. Note that once you have loaded the file, you then need to click Apply in the Properties subwindow to create the visualisation, and use the pull-down boxes near the top to select the different properties and the type of plot. To calculate new properties from the data, you can apply the Calculation filter: this allows you to create a new property from a mathematical expression based on those provided in the file. For instance, the total particle densities can be calculated using:

density_A+density_B

while the volume fraction of species A can be obtained with:

density_A/(density_A+density_B)

(Again, you will need to click Apply in the Properties subwindow to carry out the calculation.) There is also a Plot Over Line filter to produce a graph of the available properties along a given line, which can go along a particular axis.

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.

The flory_huggins.py script will automatically run a series of DL_MESO_DPD calculations to find values of \(\chi^{\text{AB}}\) for different values of \(A_{ij}^{\text{AB}}\). Full details of what it does and its defaults are given in Using DL_MESO_DPD, although for this exercise there are a few things you can do to speed up the simulations it runs (with the corresponding command line flags):

  • Decrease the dimensions of the simulation box: --L 10.0 --W 4.0
  • Increase the change in \(A_{ij}^{\text{AB}}\) between simulations: --dA 2.0
  • Run DL_MESO_DPD on multiple processor cores: --nproc 4

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 between these quantities.

Tasks

  1. Run DL_MESO_DPD in your working directory with the supplied input files. Use either traject_vtf.exe or traject_xml.exe to convert the HISTORY file data into the format to load into VMD or OVITO respectively. 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, then open Paraview and load in the averages.vtk file.
    • 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 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 (2) to calculate \(\chi^{\text{AB}}\).
  3. Either by hand or by using the flory_huggins.py Python script, 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)\).
    • Do all of the simulations produce clearly separated regions for each particle species? Which ones do not manage this? (Why might that be?)
    • Would you agree there is a proportional relationship between \(\Delta A_{ij} = \left(A_{ij}^{\text{AB}} - A_{ij}^{\text{AA}}\right)\) and \(\chi^{\text{AB}}\)? If so, what is the proportionality constant?
    • Is there a minimum value of \(\Delta A_{ij}\) that gives separation of the two species?
  4. Run some more simulations to calculate values of \(\chi^{\text{AB}}\), only this time try a different value of \(A_{ij}^{\text{AA}}\), adjusting the range of values for \(A_{ij}^{\text{AB}}\) to produce some more data points for the previous plots. Do these additional points still lie on more-or-less the same line as before?
  5. Optionally vary the system particle density from the default of 3.0 to a higher value: say either 5.0 or 6.0. (You will need the flory_huggins.py script to at least obtain the necessary CONFIG file.) 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?

Reference

[Groot1997]RD Groot and PB Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, Journal of Chemical Physics, 107, 4423-4435 (1997), doi: 10.1063/1.474784.