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 = galamost.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 = galamost.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 in GALAMOST 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 = galamost.ParticleSet(all_info, "charge")
kappa=0.2

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

# reciprocal space
enuf = galamost.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 for the input of GALAMOST. 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
sys.path.append('/opt/galamost3/lib')
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")