DPD Exercise 2: Mesophase behaviour with DPD

Introduction

The soft repulsive potentials often used in DPD simulations are one of the method’s biggest features. These allow simulation timesteps that are at least a couple of orders of magnitude larger than those available with atomistic molecular dynamics (MD). We can therefore go from randomised distributions of various components to ordered equilibrium states in relatively few timesteps.

DPD is often used to determine equilibrium structures of molecules in solution. Taking what we learned from DPD Exercise 1: DPD, hydrophobicity and parameterisation, we can join together beads with different conservative force parameters \(A_{ij}\) between these and solvent beads to produce molecules with different behaviours in different regions. Even using simple harmonic springs between these beads vastly increases the range of systems that can be modelled using DPD, including polymers, colloids and biological compounds.

The conservative force parameter indicates how much one bead is repelled by another. Many DPD simulations therefore examine amphiphiles, molecules consisting of both hydrophobic (water-hating) and hydrophilic (water-loving) sections. Surfactants are commonly amphiphilic with hydrophilic head groups and hydrophobic tails, and the structures they form depend strongly on their concentration in solution. We can therefore describe different mesophases - types of structure that the molecules can form - and construct phase diagrams using DPD simulations.

The most common mesophases are isotropic, hexagonal and lamellar. Isotropic phases tend to consist of spherical ‘blobs’ of material, hexagonal phases are made up of long tubes laid out in a hexagonal pattern, while lamellar phases consist of parallel ‘sheets’. Examples of these phases based on plots of surfaces created by hydrophobic sections of the amphiphilic molecules are shown below.

image0 image1 image2
Isotropic \((L_{1})\) phase Hexagonal \((H_{1})\) phase Lamellar \((L_{\alpha})\) phase
\(\mu_1,\mu_2,\mu_3 \approx \tfrac{1}{3}\) \(\mu_1 < 0.1 \ll \mu_2, \mu_3\) \(\mu_1 < 0.1, \mu_2 < 0.15 \ll \mu_3\)

As well as looking at trajectories, we can determine phases by calculating order parameters based on local densities of individual bead species. By assigning a Gaussian density function for each particle (on a grid) and defining a surface based on a threshold density, an isosurface normal distribution \(p(\mathbf{n})\) can be defined. The second moment of \(p(\mathbf{n})\) gives a symmetric tensor:

(1)\[\mathbf{M} = \int_{}^{}{\mathbf{nn} p\left( \mathbf{n} \right)}\text{ }d\mathbf{n}\]

and the three eigenvalues of this tensor (\(\mu_{1}\), \(\mu_{2}\) and \(\mu_{3}\)), expressed in numerical order and summed to 1, can be used to describe the phase.

Prinsen, Warren and Michels [Prinsen2002] determined guidelines for detecting isotropic, hexagonal and lamellar phases based on the values of the eigenvalues: these can be seen in the table above and can speed up phase detection considerably, making it possible to create a phase diagram using DPD simulations more quickly than real-life experiments. An example of this is given below with two qualitatively similar phase diagrams for the surfactant hexaethylene glycol monododecyl ether (often referred to as C 12 E 6): the phase diagram based on DPD simulations took about 12 hours to complete, while the other one required experiments over a few weeks.

image3 image4
Determined from DPD simulations [Seaton2011] Determined experimentally [Tiddy1980]

Aim

We are going to examine the mesophase behaviour of a simple two-bead (dimer) amphiphilic molecule in a solvent. By trying out different concentrations of the dimer in solution, we can look at the mesophases that spontaneously form even when the molecules are initially dispersed randomly in the solvent. As well as visualising the phases, we will also work out the order parameters for those phases - the eigenvalues of the second moment of isosurface normals, equation (1) - and see if we can use those to help find the boundaries between different phases.

Instructions

Copy the CONTROL and FIELD input files from the directory DEMO/DPD/AmphiphileMesophases in your copy of DL_MESO to your working directory.

These input files represent DPD simulations of a periodic cuboidal box with dimensions \(10 \times 10 \times 20\) and a total of 12,000 beads, some of which will be molecules consisting of one hydrophilic bead (species A) and one hydrophobic bead (species B) joined together with harmonic bonds. The CONTROL file is the same for all simulations and will run for 200,000 timesteps, using the first 20,000 timesteps to equilibrate the system before collecting trajectory data.

The four FIELD files (FIELD-30, FIELD-55, FIELD-75, FIELD-90) include different concentrations of amphiphilic molecules, with the number at the end of each file giving the volume fraction of molecules in the system. Note that only the numbers of beads in the solvent (species C) are specified under the SPECIES keyword, while the numbers of beads in molecules are determined from the molecule definition (using the nummols keyword). All FIELD files use the same conservative interaction parameters \(A_{ij}\) between each pair of species.

Before running each simulation, create a copy of a FIELD-* file and rename it to FIELD, e.g.

cp FIELD-30 FIELD

Afterwards, produce a traject.vtf file (or a series of traject_*.xml files) using the traject_vtf.exe (traject_xml.exe) utility.

To obtain isosurface plots of a particular species and the order parameters, you can use the isosurfaces.exe utility, which can be run as follows to generate the volume map for species B (the hydrophobic beads in the molecules):

./isosurfaces.exe -b B

There are additional command-line options to determine the grid fineness and size of Gaussian functions used to produce the volume maps, but the default values will be good enough for this exercise. As the utility runs, it will also calculate the second moments of the isosurface normal distributions based on equation (1) and display the three eigenvalues on the screen for each trajectory frame: the same numbers will also be written to a file called moment, which can be loaded and plotted using your favourite graph plotting software (e.g. Excel or Gnuplot).

Tasks

  1. Run DL_MESO_DPD in your working directory with the supplied input files for each concentration. Use either the traject_vtf.exe or traject_xml.exe utilities to obtain plottable trajectories from the HISTORY files you can load into VMD or OVITO. Can you see the various mesophases in each simulation? (I suspect it might be difficult, so don’t worry if you cannot see them!)
  2. Use the isosurfaces.exe utility to create volume maps for each simulation and calculate the order parameters. Load the generated files (density_*.vtk) into Paraview and apply the Iso Volume filter, adjusting the minimum and maximum densities until you can see distinctive shapes for the various mesophases.
    • Can you work out - from the volume plots - which type of mesophase is formed by each molecular concentration? (Hint: for a concentration of 90%, you might want to try using the isosurfaces.exe utility on the solvent - species C - rather than the hydrophobic beads.)
    • Take a look at the order parameters and compare them to the suggestions in the first table above. Do you think those parameters can easily identify the phases? How quickly does each phase form? Has each phase formed by the end of the equilibration period (the first trajectory frame)?
  3. Try modifying the molecular concentration by changing the number of molecules in the FIELD file, adjusting the number of solvent beads to ensure the total number of beads remains the same.
    • Can you find the concentrations where the various boundaries between the phases (isotropic/hexagonal, hexagonal/lamellar, lamellar/isotropic) happen to be?
    • How distinct are the phase boundaries? (Are the changes in mesophase sudden or do you get ranges of concentrations with some mixtures of phases?)
  4. Try modifying the simulation temperature in the CONTROL file and repeat the calculations: do the mesophase boundaries change with temperature?
[Prinsen2002]P Prinsen, PB Warren and MAJ Michels, Mesoscale simulations of surfactant dissolution and mesophase formation, Physical Review Letters, 89, 148302, 2002, doi: 10.1103/PhysRevLett.89.148302.
[Seaton2011]MA Seaton and JA Purton, New capabilities and software for materials modelling, p. 20-23 IN: Science Highlights 2010-2011, STFC, 2011, https://www.scd.stfc.ac.uk/Pages/8Michael_Seaton.pdf
[Tiddy1980]GJT Tiddy, Surfactant-water liquid crystal phases, Physics Reports, 57, p. 1-46, 1980, doi: 10.1016/0370-1573(80)90041-1.