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 flory_huggins.py and flory_huggins_plot.py Python scripts.

To get started, launch the flory_huggins.py script using one of the following commands:

python flory_huggins.py --dlmeso <dlmeso>
python3  flory_huggins.py --dlmeso <dlmeso>
./flory_huggins.py --dlmeso <dlmeso>

replacing <dlmeso> with the location of your DL_MESO_DPD executable, dpd.exe. With all other default settings, this script will launch 11 DL_MESO_DPD calculations with different \(A_{ij}^{\text{AB}}\) values and generate a data file (floryhuggins-rho-3.000.dat) with concentration profiles and \(\chi^{\text{AB}}\) values for each \(A_{ij}^{\text{AB}}\) value.

Full details of what this script 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 each DL_MESO_DPD calculation 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. 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)\). (The script will do this by default.)

    • 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?

  2. 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?

  3. Optionally vary the system particle density from the default of 3.0 to a higher value: say either 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?

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.