DPD Tutorial Exercise 1: Simple DPD systems

Ex1. DPD fundamentals and code hacking

Take a look at dpd.F90. You should be able to make out various sections and subroutines:

  • reading an input file (to come)

  • allocations of arrays for positions, velocities and forces of particles

  • initial simulation setup (placing particles and assigning velocities)

  • opening output files for statistics and trajectory data

  • a main loop through the timesteps:

    • applying stage 1 of Velocity Verlet integration of forces

    • calculating forces on pairs of particles within the cutoff distance

    • applying stage 2 of Velocity Verlet

    • writing to output files and to the screen at specified intervals

  • calculations and printing of final properties

  • closing of files and the program itself

Try the following command to compile the code:

[ ]:

!gfortran -o DPD dpd.F90

You should find that the code does not compile and throws up error messages for lines 603 and 609: this is expected behaviour, so please do not panic! These two lines are the ones meant to calculate the dissipative and random forces for a given pair of particles, but at the moment they are incomplete. Your main task for this exercise is to complete these lines and try out the code for a single-component DPD simulation.

Some hints on completing those lines:

  • The dissipative force parameter \(\gamma\) is available as gamma in the code.

  • The variable sigma is equal to \(\sqrt{\frac{2\gamma k_B T}{\Delta t}} = \frac{\sigma}{\sqrt{\Delta t}}\), i.e. the random force parameter divided by the square root of the timestep size. (This reduces the number of divisions the code has to make.)

  • The simplest option for the screening function for random forces \(w^R\) would be the equivalent of that for conservative forces, given as wrr, and a Gaussian random number \(\zeta\) is already calculated for you, available as zeta.

  • The screening function for dissipative forces must equal \(\left(w^R\right)^2\), and the dot product of the vector between the particles and their relative velocity, \(\textbf{r}_{ij} \cdot \textbf{v}_{ij}\), is available as rdv.

  • The forces fdr and frr need to be divided by \(r_{ij}\) (rrr) as these will then be multiplied by the actual vector between the particles: see lines 545-550. Note that you may have to divide fdr by rrr twice, and \(r_{ij}^2\) is available as rsq.

Once you are done, try compling the code again:

[ ]:

!gfortran -o DPD dpd.F90

Take a quick look at the input file we are going to use with this code: INPUT. Then run the code using the following command to ‘pipe in’ the input file:

[ ]:

!./DPD < INPUT

If the code has run as expected, you should see some checkpoint data on the screen in four columns: current time in DPD units, total system energy, temperature and pressure. You will also end up with three output files:

  • STATS - a statistics file with time, system energy, temperature and pressure for the entire simulation run in tabulated (fixed-width) columns

  • XYZ - the final configuration of the system in XYZ format

  • TRAJECT - the system trajectory: a sequence of configurations taken at intervals during the simulation in XYZ format

You can take a look at the contents of STATS by using the commands below to invoke a simple plotting script or by opening the file in plotting software (e.g. Gnuplot, Xmgrace) or a spreadsheet program. XYZ and TRAJECT can be opened and visualised in VMD, selecting XYZ as the file format in each case.

[ ]:

from plotstats import *
plotStats('STATS', 1, 2, 'Time', 'Energy', 'Energy vs. time')
plotStats('STATS', 1, 3, 'Time', 'Temperature', 'Temperature vs. time')
plotStats('STATS', 1, 4, 'Time', 'Pressure', 'Pressure vs. time')

Taking a look at these files:

  • Is the simulation stable? Does the temperature correspond to the value in the input file?

  • What happens to the particles when they get close together?

  • Try re-running the program. You are unlikely to get exactly the same results due to how the random number generator is seeded in the code, but how close are the energies, temperatures and pressures between runs?

  • What happens to the temperature if you halve or double the timestep size?

In case you struggle to correctly modify the code, we have a correct version of it available to check: dpd_solution.F90.

So far, we have used a brute force approach for finding particle pairs within the cutoff distance, i.e. two nested loops and an IF statement to check the distances. A more sophisticated linked-cell list approach is also available, which divides the simulation box into cubes (‘cells’) with sides of at least the cutoff distance, assigns each particle into a cell and searches for particle pairs within each cell and its nearest neighbours but no further!

We have programmed the linked-cell list approach into our code and used preprocessor directives, which means we can use a command-line option with the Fortran compiler to invoke different parts of the code. To compile our DPD code with this option, try the following command:

[ ]:

!gfortran -o DPD-FAST -DFAST dpd.F90

and then run the code with the same input file using the following command:

[ ]:

!./DPD-FAST < INPUT

It should generate the same files as before, which you can open in the same ways as mentioned above.

  • For the same system size (number of particles and volume), do you notice any difference in reported calculation times?

  • Try increasing the number of particles and volume by the same factor (keeping the density constant) in the input file, and run the simulation with both versions of the code. How do the calculation times change for each version of the code compared with the original system size?

  • Can you work out how well each code version scales with the number of particles?

The code is not especially exciting, but it can allow you to explore some properties of a DPD fluid. If you want to explore these, try the following optional tasks:

  • Vary the temperature and particle density to see if you can obtain a stable crystalline structure.

  • Use the TRAJECT file in VMD to calculate diffusion constants for higher temperature systems.

  • Try varying the ‘energy parameter’ (conservative force parameter \(A_{ij}\)) and measure the fluid pressure at different temperatures and densities.

  • See if you can find values of the ‘energy parameter’ that will represent water at various coarse-graining levels (number of molecules per particle).

The next exercise is available in this notebook.