Dissipative particle dynamics¶

DPD force¶

Description:

The DPD force consists of pair‐wise conservative, dissipative and random terms.

\begin{eqnarray*} \vec{F}_{ij}^{C}&=&\alpha\left(1-\frac{r_{ij}}{r_{cut}}\right)\vec{e}_{ij} \\ \vec{F}_{ij}^{D}&=&-\gamma\omega^{D}(r_{ij})(\vec{e}_{ij} \cdot \vec{v}_{ij} )\vec{e}_{ij} \\ \vec{F}_{ij}^{R}&=&T\sigma\omega^{R}(r_{ij})\xi_{ij}\vec{e}_{ij} \\ \end{eqnarray*}

\(\gamma=\sigma^{2}/2k_{B}T\)

\(\omega^{D}(r_{ij})=[\omega^{R}(r_{ij})]^2=(1-r_{ij}/r_{cut})^2\)

\(\xi_{ij}\) - a random number with zero mean and unit variance

\(T\) - temperature - optional: defaults to 1.0

\(r_{cut}\) - r_cut (in distance units) - optional: defaults to 1.0

The following coefficients must be set per unique pair of particle types:

\(\alpha\) - alpha (in energy units)

\(\sigma\) - sigma (unitless)

class DPDForce(all_info, nlist, r_cut, temperature, rand_num)¶

The constructor of a DPD interaction object.

Parameters:

all_info (AllInfo) – The system information.
nlist (NeighborList) – The neighbor list.
r_cut (float) – The cut-off radius.
temperature (float) – The temperature.
rand_num (int) – The seed of random number generator.

class DPDForce(all_info, nlist, r_cut, rand_num)¶

The constructor of a DPD interaction object. The default temperature is 1.0.

Parameters:

all_info (AllInfo) – The system information.
nlist (NeighborList) – The neighbor list.
r_cut (float) – The cut-off radius.
rand_num (int) – The seed of random number generator.

setParams(string typei, string typej, float alpha, float sigma)¶: specifies the DPD interaction parameters per unique pair of particle types.

setT(float T)¶: specifies the temperature with a constant value.

setT(Variant vT)¶: specifies the temperature with a varying value by time step.

setDPDVV()¶: calls the function to enable DPDVV method (the default is GWVV).

Example:

dpd = gala.DPDForce(all_info, neighbor_list, 1.0, 12345)
dpd.setParams('A', 'A', 25.0, 3.0)
app.add(dpd)

GWVV integration¶

Description:

Integration algorithm.

\begin{eqnarray*} &v_i^0&\leftarrow v_i+ \lambda\frac{1}{m}(F_i^c \Delta t + F_i^d \Delta t + F_i^r \sqrt{\Delta t})\\ &v_i&\leftarrow v_i+ \frac{1}{2}\frac{1}{m}(F_i^c \Delta t + F_i^d \Delta t + F_i^r \sqrt{\Delta t})\\ &r_i&\leftarrow r_i+ v_i \Delta t\\ &&Calculate \quad F_i^c\{r_j\}, F_i^d\{r_j, v_j^0\}, F_i^r\{r_j\}\\ &v_i&\leftarrow v_i+ \frac{1}{2}\frac{1}{m}(F_i^c \Delta t + F_i^d \Delta t + F_i^r \sqrt{\Delta t}) \end{eqnarray*}

\(\lambda\) - lambda (unitless) - optional: defaults to 0.65

class DPDGWVV(AllInfo all_info, ParticleSet group)¶

The constructor of a GWVV NVT thermostat for a group of DPD particles.

Parameters:

all_info (AllInfo) – The system information.
group (ParticleSet) – The group of particles.

setLambda(float lambda)¶: specifies lambda.

Example:

gwvv = gala.DPDGWVV(all_info, group)
app.add(gwvv)

Coulomb interaction in DPD¶

Description:

In order to remove the divergency at \(r=0\), a Slater-type charge density is used to describe the charged DPD particles. Thereby, Ewald for DPD (short-range) (DPDEwaldForce) method can be employed to calculated the short-range part of Ewald summation. The long-range part of Ewald summation can be calculated by PPPM (long-range) (PPPMForce) or ENUF (long-range) (ENUFForce). And the ENUF (long-range) (ENUFForce) is suggested.

Example:
group_charge = gala.ParticleSet(all_info, "charge")
kappa=0.2

# real space
ewald = gala.DPDEwaldForce(all_info, neighbor_list, group_charge, 3.64)#(,,rcut)
ewald.setParams(kappa)
app.add(ewald)

# reciprocal space
enuf = gala.ENUFForce(all_info, neighbor_list, group_charge)
enuf.setParams(kappa, 2, 2, 20, 20, 20)
app.add(enuf)

Reduced charges:

The charges should be converted into the ones in reduced units according to Charge units. Typically, the fundamental length and energy are \(\sigma=0.646\text{ }nm\) and \(\epsilon=k_{B}T\) with \(T=300\text{ }K\), respectively, in DPD. The reduced charges are \(q^{*}=z\sqrt{f/(\sigma k_{B}T \epsilon_r)}\). The \(z\) is the valence of ion.

Here is a molgen script for polyelectrolyte.

Example:
#!/usr/bin/python
import sys

import molgen
import math

er=78.0
kBT=300.0*8.314/1000.0
r=0.646
gama=138.935
dpdcharge=math.sqrt(gama/(er*kBT*r))

mol1=molgen.Molecule(50)
mol1.setParticleTypes("P*50")
topo="0-1"
for i in range(1,50-1):
  c=","+str(i)+"-"+str(i+1)
  topo+=c
mol1.setTopology(topo)
mol1.setBondLength(0.7)
mol1.setMass(1.0)

mol2=molgen.Molecule(1)
mol2.setParticleTypes("C")
mol2.setMass(1.0)
mol2.setCharge(dpdcharge)

mol3=molgen.Molecule(1)
mol3.setParticleTypes("A")
mol3.setMass(1.0)
mol3.setCharge(-dpdcharge)

mol4=molgen.Molecule(1)
mol4.setParticleTypes("W")
mol4.setMass(1.0)

gen=molgen.Generators(15,15,15)
gen.addMolecule(mol1,1)
gen.addMolecule(mol2,75)
gen.addMolecule(mol3,75)
gen.addMolecule(mol4,9925)

gen.outPutXml("ps0")