| Property Prediction |
Sections in this chapter
parameters as input to DPD simulations
In the DPD methodology, the fundamental particles are "beads" that represent small regions of fluid material rather than the atoms and molecules familiar from MD simulations. All degrees of freedom smaller than a bead radius are assumed to have been integrated out leaving only coarse-grained interactions between beads. There are three types of force between pairs of beads, each of which conserves both bead number and linear momentum: an harmonic conservative interaction, a dissipative force representing the viscous drag between moving beads (i.e., fluid elements), and a random force to maintain energy input into the system in opposition to the dissipation. All forces are short-ranged with a fixed cut-off radius. By a suitable choice of the relative magnitudes of these forces, a system can be shown to evolve to a steady-state that corresponds to the Gibbs Canonical ensemble. Integration of the equations of motion for the beads generates a trajectory through the system's phase space from which all thermodynamic observables (e.g., density fields, order parameters, correlation functions, stress tensor, etc) may be constructed from suitable averages. An immense advantage over conventional Molecular Dynamics and Brownian Dynamics simulations is that all forces are "soft" allowing the use of a much larger time-step and correspondingly shorter simulation times.
Introduction
Using DPD
| At any time when running DPD, online help is available by clicking the right mouse button when the cursor is over a panel or a command. |
| Start with a new session of Cerius2. |
| From the Cerius2·Visualizer card deck menu go to the MESOSCALE card deck and choose the DPD card. |
This reveals a menu consisting of Run, Build, Analysis and Job Control and Reset.
| First, make sure that you have a molecular ensemble defined for the DPD simulation. From the Build submenu choose Beads. |
| Check that three beads, named A, B and W are defined there. |
| Next, from the Build submenu, choose Molecules. |
| Check that one "triblock" molecule consisting of the two beads A and B is set (i.e., A 3 B 9 A 3). |
Next you conclude the definition of the molecular ensemble.
| Close all open panels. |
| Go to the DPD card and select the Run item. |
There you find the main parameters for a particular run.
This is all you need to do for now. The control panels accessible from the Run panel provide further more detailed control which does not concern us for now.
| Go to the DPD card and select the Job Control item. |
| Finally, go to the DPD Run panel and press the RUN button. |
You should get the following message in the textport:
DPD job FirstTest has been started.On the Job Control panel the job is listed in the Job Status box.
You can monitor the progress of the job by clicking the Monitor status file button. This opens a window displaying the bottom of the status file.
| When the number of steps you had selected has been reached, press the UPDATE Job status button. |
The Status indication in the Job Control Box changes to complete, if the job has finished.
| For that purpose, choose System from the Analysis submenu on the DPD card. |
| The Morphology entry on the Analysis submenu shows two pullright options: Isodensities and Profiles. Choose Profiles first. |
You will see the name of the selected system, the time step and bead name.
| Press Create New Profile. |
| You may not see much evidence of phase separation for this short run, but you can enhance the contrast by opening the More Editing Options... control panel and pressing Reset Map to Full Range. |
The density boxes show the lowest and highest concentration of A in the slice.
| Move the slice by using the sliders. |
You can do further analysis by looking at the isodensity surfaces.
| Select the Analysis/ Morphology/Isodensities item. Set the isodensity of A to 0.4. Make sure that you are in an empty model space, and then click the Create Isodensity Display button. |
The isodensity value must fall between the minimum and maximum values you found in the Profiles analysis of that bead.
A graphical user interface module to the DPD simulation has been provided within the Cerius2 molecular modeling environment. This interface handles the input specification (i.e., the definition of the beads and polymers and their interactions, setting of the simulation run parameters, control of the numerical algorithm parameters), the job control (choice of machine for run, and output paths, job monitoring), the output file handling from the (possibly remote) hosts, and the analysis of the simulation output (3-D density fields and scalar thermodynamic functions). In addition, Cerius2 provides access to a range of molecular modeling tools that can furnish the Flory-Huggins parameters required for calculating the repulsion parameters between beads in DPD. The relationship between these two sets of parameters is described in detail in the preceding Theory section.
General methodology
|
The first action in many cases will be to build a molecular ensemble. This is done from the Build sub-menu. The Beads, Molecules, Interactions, Dissipations, and Constraints menu items are defined here (see Figure 16).
The Run panel (Figure 17) then provides for editing of the parameters governing the actual simulation, such as the length of run, time step-size, density, shear rate (if required), output frequency and the name of the run.
In the Analysis panel a particular system for which output is available can be selected and its output files inspected (Figure 19), output functions plotted (Figure 20), and the three-dimensional density profiles of the beads be analysed graphically by slicing (Figure 18) and isodensity surfacing (Figure 23).
The Job Control section (Figure 19) allows the host for the run, and its associated parameters to be set, the job to be monitored, and output to be retrieved from remote hosts for analysis. It also provides for two modes of simulation: Interactive mode, in which the DPD simulation continuously updates the Cerius2 model window with snapshots of the system, and Background mode, in which the simulation runs independently of Cerius2. Advanced options here allow the user to modify the paths to the background executable and the spawner program required for licensing purposes, as well as the username and password for runs on remote hosts.
Building a molecular ensemble
A DPD simulation consists of a system of beads which may be connected together to form polymers. The number and type of beads and their connectivity into polymers together with the forces between them are specified in the sub-panels of the Build menu item.
|
The Build menu item on the DPD card provides access to the Beads, Molecules, Interactions, Dissipations and Constraints panels which allow all of the above parameters to be set (Figure 16). In particular, the chain architecture is entered in the form of a \Qtree string' by which even complex branched structures can be defined as input to DPD. Comprehensive on-line help is available for this task. The spring constant for the Hookean force law between beads connected into a polymer is set in the Molecules panel, as are the concentrations of each molecule type (including single-bead species such as water). The interaction parameters between all pairs of beads, and between beads and the wall if it is selected in the Constraints panel, are entered here as well as the dissipative forces for all bead pairs.
Specifying the Run parameters for a DPD simulation
Having defined a molecular ensemble, a number of parameters need to be set to define the particular simulation to be carried out. These parameters are specified in the DPD Run panel (Figure 17). Each run is identified by a seed name, given as the file prefix to all input and output files generated. The model system is further defined by specifying its temperature, and the size of the simulation box (with periodic boundary conditions), in terms of the side length of its cubic cells and the number of cells in each cartesian dimension. Furthermore, the step size for the numerical integration and the length of the run, either in terms of its total time, or in terms of the number of steps, must be set. A title for the simulation may be added, to include comments about the specific run. The date and time are added to this automatically.
The output level of the run can be set from the Output... control panel. In particular, the period, in units of time steps, with which the output is sampled and the status information, density files, and restart files are written to disk is set here. The amount of disk space required for frequent output of three-dimensional fields, such as contained in density, and restart files, may be very large.
|
A simulation is started by pressing the RUN button.
|
Setting the host machine, job monitoring and output handling
Opening the Job Control panel (Figure 19) from the DPD menu shows the current host used for running DPD and the mode of execution (interactive or background). The Options... control panel allows you to set the paths to the background executable code and the spawner program required to run the licensed DPD code. Access information to the remote host (user name and password) can also be set here.
|
Further on the Job Control panel, you can set the base directory on the remote host. The Monitor box shows previous and current jobs and related information such as process id and status.
Analysing the result: Selecting a system, inspecting the files, plotting the thermodynamic functions, and exploring the morphology
Having completed a simulation and transferred the output files to the local host, you are ready for analysis of the results. This functionality is accessed via the Analysis submenu on the main DPD card.
|
As a next step, you can inspect the time evolution of the diffusion coefficients for all beads in the system by selecting the PLOT/DIFFUSION menu item (Figure 21). Alternatively, you can display the polymer endpoint distribution or the bond-length distribution for each polymer type.
|
Finally, and usually most importantly, the morphology of the molecular ensemble can be investigated and analysed. Two methods are provided for this.
|
On the DPD Profiles panel, you may change the position and direction of the selected slice by means of the sliders. More Editing Options... shows the way in which the density is mapped to color, and lets you edit this mapping. The Profile grid scaling is a factor by which the actual grid is multiplied to generate the triangular mesh for the graphics. A higher number gives a smoother appearance but takes longer to process.
|
You can open this panel from the Morphology menu item of the Analysis menu. For each bead type you can create an isodensity surface at a number of time steps and for up to four different values of the concentration at once. Typical analyses involve displays of (a) the different beads, each at a concentration close to its maximum (to be determined from the slice display), (b) one bead type at a number of time steps, and (c) one particular bead at different isodensities. Each surface in the current model can be edited in terms of its visibility, transparency and color.
Theory
Introduction
The most accurate method of calculating the dynamical behavior of an atomistic system is to integrate the equations of motion of all the atoms in the system. This is the basis of the Molecular Dynamics simulation technique. Its major drawback is that it often provides far more detail of small-scale fluctuational motion of atoms than is necessary for an understanding of many physical processes. It also requires computer processor speeds and memory capacities that currently limit its applicability to a few nanoseconds of molecular motion. This is inadequate for many chemical processes that occur on the microsecond (or longer) time-scales. Some years ago, Hoogerbrugge and Koelman (1992) introduced a new simulation technique, derived from Molecular Dynamics simulations and Lattice Gas Automata, that effectively opens up the mesoscopic length and time regimes in complex fluids to simulation. It achieves this goal by keeping the principle of integrating the equations of motion for a system but integrating out the smallest spatial degrees of freedom first. The fast motion of the atoms in a system is averaged over and the remaining structure is represented by a set of "beads", of given mass and size, that interact via soft potentials with other beads. A bead represents a small region of fluid matter and its motion is assumed to be governed by Newton's laws in which the total force on a bead is a sum of its direct interactions and dissipative and random forces between it and other beads. The dynamical behavior of the system is followed along a trajectory through its phase space by integrating its equations of motion. The equilibrium properties are calculated by performing suitable averages along this trajectory. In this section we describe in detail the simulation technique and how it can be related to the Flory-Huggins theory of polymeric fluids. The following description of the theory is based on a paper by Groot and Warren (1997).
Equations of motion for a DPD system
We consider a set of beads interacting by specified forces and whose dynamical evolution is governed by Newton's laws
Eq. 21
where ri, vi and fi and are the position vector, velocity and total force on the ith bead. For simplicity, we set all bead masses mi equal to unity. Each bead is subject to three forces from its neighbors: a conservative interaction which is linear in the bead-bead separation; a dissipative force proportional to the relative velocity of two beads and a random force between a bead and each of its neighbors. The total force on a bead is 
Eq. 22
where the sum is over all beads within a distance rc of the ith bead. This short-range cut-off makes the interactions local. From now on we set rc equal to unity so that all lengths are measured relative to the bead radius. The conservative force is a repulsive central force with a maximum magnitude aij
Eq. 23
where rij is the magnitude of the bead-bead vector rij , and r-hatij is the unit vector joining beads i and j.
Eq. 24
where 
D(rij) is a short-range weight function. Because of the form chosen for the dissipative force it conserves the total momentum of each pair of particles, and hence also of the system. The random force also acts between all pairs of beads subject to a similar short-range cut-off, with a possibly different function R(rij), and acts so as to pump energy into the system
Eq. 25
where 
ij(t) is a delta-correlated stochastic variable with zero mean
Eq. 26
Note the unusual property of the random force that it is pairwise central (in contrast to Brownian Dynamics in which white noise is added to the overdamped equation of motion for each particle independently) and hence conserves total linear momentum even while adding energy to the system. At this stage, there are two unknown functions, 
D(rij) and R(rij), and two unknown constants 
and 
. Espagnol and Warren (1995) have shown that in order for the steady-state solution to the equations of motion to be the Gibbs ensemble, and for the fluctuation-dissipation theorem to be satisfied, only one of the two weight functions D(rij), R(rij) can be chosen arbitrarily while this choice fixes the other; and the two multiplicative constants , are related by the temperature. That is, we must impose the relations
Eq. 27
where T is the absolute temperature and kB is Boltzmann's constant. For simplicity we choose
Eq. 28
At this point a word about the role of units in DPD simulations is in order.
Eq. 29
The results of a single DPD simulation may thus correspond to many physical systems depending on the values chosen for the bead mass and radius, the temperature and the interaction and dissipation parameters.
Integration scheme in DPD
A modified velocity-Verlet algorithm has been used to integrate the equations of motion. In this scheme, the current values of position, velocity and the force on a particle are used to calculate the position and velocity of the particle at the next time-step; the new position and velocity are then used to calculate the new force, and this then corrects the velocity.
Eq. 30
The factor 
is set to 1/2 in the DPD simulation code, although Groot and Warren (1997) discuss the effect of alternative values for this factor in taking the stochastic nature of the force into account when calculating the order of the integration scheme.
Choosing the dissipation and random noise magnitudes
It was stated above that the magnitudes of the dissipative and random forces are connected by the fluctuation-dissipation theorem. This still leaves one of them as a free parameter. We take this to be the magnitude of the dissipative forces in the system. Groot and Warren (1997) found that there was no observable difference between the effects of Gaussian noise and uniform noise and so the simpler uniform noise has been used. When the noise magnitude was larger than approximately = 8, they found that the integration scheme was unstable. It was found that a noise amplitude of = 3 gave reasonably fast relaxation for temperatures between KBT = 1 and KBT = 10. When choosing values for the dissipation parameters, care must be taken that the chosen temperature and dissipation parameters do not result in noise that is much larger than = 3 or the simulation will be unreliable.
Choosing the repulsion parameters
Once the mass and radius of the beads, the architecture of the polymers in the system, and the physical constraints of temperature, dissipation magnitude and system size have been fixed, there remains only one parameter to characterize the interactions between the beads: the repulsion parameters aij. These are all that appear in the DPD methodology to represent the complex forces that act between the atoms and molecules in a physical system.
Eq. 31
where p is the pressure, 
is the density and 
= 0.101 ± 0.001. In order to fix the magnitude of the repulsion parameters we compare the dimensionless compressibility of the DPD fluid with that of water using the equation
Eq. 32
For water at room temperature (300 K) the dimensionless compressibility has the value 
-1 = 15.9835. By differentiating the equation of state and comparing with the compressibility of water we find a/kBt 
75¾. In principle, we could choose the density freely, but to minimise the number of interactions between beads in a simulation we take the smallest value at which this relationship holds reasonably well: Groot and Warren have shown that 3 is sufficient. In order that a DPD fluid has the compressibility of water, we need a repulsion parameter of a 25kBT at 3; for other densities we use a 75kBT/.
Mapping the interactions onto Flory-Huggins theory
Having established that a DPD simulation can represent the volume fluctuations of a simple fluid, we would like next to be able to model the interactions between unlike fluids such as copolymer melts or surfactant-oil-water mixtures. We need to find a relation between the interactions between DPD beads and the theory of polymer mixtures so that the behaviour observed in a simulation can be mapped onto the physical phase diagram of real fluids. One way of doing this is to compare the free energy of a DPD fluid with the mean-field theory of polymer mixtures, Flory-Huggins theory. We give a brief summary of this theory.
A and B the volume fractions of each polymer. The free energy of the mixture in excess of the ideal gas contribution is
Eq. 33
where the filled lattice condition means 
A + B = 1. The parameter controls the interaction between the polymers: if positive, the two species prefer to phase separate, whereas if it is negative they prefer to mix. The equilibrium state of the mixture is found by minimising the free energy which leads to an implicit relation for the density
Eq. 34
As 
increases from a negative value the system starts to phase separate. The critical value of at which phase separation first occurs may be found by minimising the free energy and noting where the spinodal points coincide. As in van der Waal's theory of gases, this occurs when the first and second derivatives of the free energy are zero. We find
Eq. 35
The free energy density of a single-component DPD fluid whose equation of state, as given above, is quadratic in the density is
Eq. 36
and for a two-component fluid we expect
Eq. 37
Setting aAA = aBB and assuming that the total density is approximately constant, we find (up to constant terms)
Eq. 38
where we have written x = 
A/(A + B), and identified
Eq. 39
Comparing the DPD free energy with the Flory-Huggins equivalent suggests that the two theories correspond if the 
parameter is proportional to the repulsion parameters as specified by this equation. Groot and Warren (1997) find that the excess pressure of a binary mixture of DPD monomers and polymers is proportional to the term x(1-x) for repulsion parameters aAA = aBB = 15, 25, , but that the constant of proportionality is not linear in the repulsion parameters (aAB = aAA). However, for large enough values of the parameter we expect the mean field theory to be reasonably accurate, and we can use the above relation between the density of the DPD fluid and the equivalent parameter to relate the two theories. The average density is measured in two homogenous regions, distant from any interfaces, in a DPD simulation of a binary mixture. The quantity NA is obtained from the density via the logarithmic expression above (noting that NA = 1 for monomeric systems), and is plotted against the repulsion parameter difference, (aAB = aAA) to test the correspondence between the Flory-Huggins theory and the simulation. Groot and Warren (1997) find
Eq. 40
for densities of 
= 3.5 respectively. Notice that although the coefficient is not linear in the density as expected, such a relation can be measured from a set of simulations at a fixed density and used to relate the repulsion parameters in the DPD fluid to the parameter in the Flory-Huggins theory. are chosen to be
Eq. 41
where the second term is absent for interactions between like beads. The 
parameter for polymers is found from simulations by varying the repulsion parameters and the polymer length and plotting NkBT/(aAB - aAA) against N. Groot and Warren (1997) find
Eq. 42
This relation enables the properties of real polymers, with lengths of the order of ten thousand monomers, to be simulated using polymers with only 10 or so DPD beads per polymer, if the 
parameter from the Flory-Huggins theory of the polymers is increased such that the product N remains constant. The input to the DPD simulations may be obtained from MD simulations of realistic models of the polymers that generate only the parameters. We next indicate how these may be obtained from microscopic simulations and other techniques within Cerius2.
Calculation of Flory-Huggins
There are several methods of calculating the parameters as input to DPD simulations parameters for polymers of known structure. Three of these are available within the Cerius2 software suite. In order of increasing complexity they are:
The second method relies on calculating the free energy of mixing of one species of polymer in another from their pair contact energies. Flory-Huggins theory assumes that the
parameter for a mixture of A and B polymers represents the repulsive energy of an AB pair averaged over all possible configurations: the pair is considered to feel only the mean field of all the other polymers and to adopt that configuration that minimises the mean field free energy. Hence, it is natural to try and calculate the pair contact energy by sampling over the most probable conformations of a pair of polymers in contact using Monte Carlo simulations. By suitably choosing the set of trial conformations, the average interaction energy can be calculated and the parameter extracted.The final method involves using Molecular Mechanics simulations to find the cohesive energy of the two species of polymer in their pure phase and when mixed together. Comparing these with the gas phase energy of each species allows us to calculate the energy of mixing of the polymers and extract the
parameter. This method is more accurate than the other two, including as it does hydrogen bonding effects and chain structure in atomistic detail, but takes considerably longer to calculate.Further details on these methods can be found in the Cerius2 documentation for the relevant products.
Case, F. H. and D. J. Honeycutt. TRIPS. 2 259 (1994).
Espagnol, P. and P Warren. Statistical mechanics of Dissipative Particle Dynamics (1995).
Europhys. Lett. 30 191-196.
Fan. C. F., B. D. Olafson, M. Blanco, and S. L. Hsu. Macromolecules. 25 366 (1987)..
Groot R. D. and P. B. Warren. "Dissipative Particle Dynamics: Bridging the gap between atomistic and mesoscopic simulation" J. Chem. Phys. 107 4423-4435 (1997).
Hoogerbrugge, P. J. and J. M. V. A. Koelman. "Simulating microscopic hydrodynamic phenomena with Dissipative Particle Dynamics" Europhys. Lett. 19 155-160 (1992).