DPD Tutorial Exercise 2: Applying DPD to molecular systems¶
Ex2. Lipid bilayers, micelles and vesicles¶
Following on from amphiphilic dimers in the previous exercise, we can construct longer, more realistic amphiphilic molecules - often with smaller hydrophilic head groups and longer hydrophobic tails - such as surfactants and lipids (naturally occurring hydrophobic or amphiphilic molecules, e.g. fats, waxes, sterols). Larger molecules can incorporate more complex topologies, such as interactions that control the angles between adjacent bonds.
In this exercise, we will use a DPD model for a lipid molecule consisting of a hydrophilic bead and six hydrophobic beads joined together with harmonic bonds, along with a cosine bond angle potential:
to conctrol the angle between pairs of bonds. The angles will have an effect on the structures that the lipids form, including:
micelles - spherical structures with hydrophobic tails pointing inwards
bilayers - two planar layers of lipids with hydrophobic tails pointing towards each other
vesicles (liposomes) - spherical kinds of bilayer with hydrophilic head groups and solvent inside
Since the required simulations are fairly large, we are going to use SCARF to run DL_MESO_DPD. We have already compiled the parallel version of DL_MESO_DPD and one of the utilities (traject_vtf.exe
) used to convert HISTORY
files to VTF files for VMD to run on SCARF. We have also uploaded input files (CONTROL and FIELD) and created a Slurm job script to run our lipid simulation. Note that most of the following require you to
use a terminal window, and that you have already setup access to SCARF.
Log on to SCARF:
ssh scarf
and create your work directory if it does not yet exist (substituting tmpXY
with your assigned username):
mkdir /work4/training/ccp5/tmpXY
Go into that directory (cd /work4/training/ccp5/tmpXY
) and then copy over the input files and job script into it from a shared directory:
cp ../meso/CONTROL .
cp ../meso/FIELD .
cp ../meso/dpdjob .
We want to run this simulation four times, so create four sub-folders and copy all three files into each of them, e.g.
mkdir LIPID1 LIPID2 LIPID3 LIPID4
cp CONTROL FIELD dpdjob LIPID1/
cp CONTROL FIELD dpdjob LIPID2/
cp CONTROL FIELD dpdjob LIPID3/
cp CONTROL FIELD dpdjob LIPID4/
Simulation 1¶
Go into the first folder and launch the simulation using the job script:
cd LIPID1
sbatch dpdjob
This will put your simulation into the job queue and, when there are sufficient processor cores available, run DL_MESO_DPD on 32 cores for up to an hour before running the traject_vtf.exe
utility to generate a traject.vtf
file from the HISTORY
file. (We have recently checked that the simulation should take around 10-15 minutes to run.) Note that while the next few steps that involve looking at the results will require you to wait until this simulation has finished, the other
simulation setup steps will not, so feel free to skip ahead and come back here later after launching the other three simulations. (If the job queue is reasonably empty, all four calculations might start and finish at around the same time.)
Once the simulation has finished running and the traject.vtf
file has been created, type the following command on SCARF inside the directory to find its full path:
pwd
and then open a new terminal instance on here, navigate to the Day_9Meso
directory and then invoke the following commands to copy the traject.vtf
, CORREL
and OUTPUT
files into the (intentionally initially empty) DPD2Ex2/LIPID1
directory, replacing /path/to
with whatever the above pwd
command gave you:
scp scarf:/path/to/traject.vtf DPD2Ex2/LIPID1
scp scarf:/path/to/CORREL DPD2Ex2/LIPID1
scp scarf:/path/to/OUTPUT DPD2Ex2/LIPID1
Take a look at the traject.vtf
file in VMD, the OUTPUT
file and the CORREL
file using the script below:
[ ]:
from plotcorrel import *
plotCorrel('DPD2Ex2/LIPID1/CORREL', 'time', 'pe-total', 'Time', 'Potential energy per particle', 'Potential energy per particle vs. time')
plotCorrel('DPD2Ex2/LIPID1/CORREL', 'time', 'temperature', 'Time', 'Temperature', 'Temperature vs. time')
plotCorrel('DPD2Ex2/LIPID1/CORREL', 'time', 'bndlen-av', 'Time', 'Mean bond length', 'Mean bond length vs. time')
plotCorrel('DPD2Ex2/LIPID1/CORREL', 'time', 'angle-av', 'Time', 'Mean bond angle (degrees)', 'Mean bond angle vs. time')
What structure has formed during this simulation?
Does the formation of this structure coincide with the lowest potential energy levels for the simulation?
What averaged bond angles do you end up getting?
Simulation 2¶
For our second simulation, we will try and ‘switch off’ the angle potentials. Go into your second sub-folder on SCARF (LIPID2
), open the FIELD
file in a text editor and find the following section defining the angle potentials for the molecule:
angles 5
cos 1 2 3 20.0 0.0 1.0
cos 2 3 4 20.0 0.0 1.0
cos 3 4 5 20.0 0.0 1.0
cos 4 5 6 20.0 0.0 1.0
cos 5 6 7 20.0 0.0 1.0
The last three columns in these lines define \(A\), \(\theta_0\) and \(m\) in our cosine angle potential. For this simulation, change the values of \(A\) from \(20.0\) to \(0.0\), then save the FIELD
file and launch the calculation. Copy over the output files into the DPD2Ex2/LIPID2
folder and take a look at them as you did for the first simulation.
What structures have the lipids now formed?
What bond angles do you now get in this simulation?
Simulations 3 and 4¶
For our third and fourth simulations, we can try modifying the concentration of molecules to see how this affects the structures that the lipids form.
Go into your third sub-folder on SCARF (LIPID3
), open the FIELD
file and reduce the number of HC6
molecules (the value after nummols
) from 2000 to a lower number: say between 800 and 1500. While you have the same file open, increase the number of solvent (W
) beads to compensate for the reduced number of molecules and to keep the total number of beads constant at 49152. (The volume of the box is \(16384.0\) and we want a bead density of \(\rho = 3\).)
Make the same changes to the FIELD
file in the fourth sub-folder but also change the value of \(A\) for the cosine bond potentials to \(0.0\) as you did for the second simulation.
Launch both calculations and copy over the output files into the DPD2Ex2/LIPID3
and DPD2Ex2/LIPID4
folders respectively so you can take a look at them.
What effect does reducing the molecular concentration have on the lipid structures, compared with the first and second simulations?
Do you observe any noticeable differences between the third and fourth simulations, given the fourth has switched off the angle potentials?