Problems_DPD_parameters

I am initiating a studying case of micellization of a nonionic surfactant in water. The simulations were terminated in a few steps due to the divergent energy, and I figured it could be my misunderstandings for some parameters introduced in the LAMMPS manual. What I did is addressed as follows. Input file and the data file are described in the bottom.

The system contains two molecule types, surfactant and water. A surfactant is modeled by 6 beads (1-1-2-2-2-2) containing two bead types 1 and 2, and water is modeled as bead type 3. Bead size (Rc) is unified as 7.1 nm, and all length units in the system are reduced by this value (box length, bead coordinates, bond length, etc.) Bead mass is unified as 1. The simulation box is 30x30x30 (Rc^3). Surfactant bond types (two body for connection and three body for rigidity) are assigned using harmonic bonds. Number of beads is 81000, projecting to the density equals to 3 in the box. Gamma (force/velocity unit) in DPD is set to 1.5. Initial coordinates are taken from the final configuration of my previous calculations under same conditions. After performing parallel simulations using 8 cores, I had error message “Bond atoms %d %d missing on proc %d at step %ld.” I think that I must set the wrong pairwise parameters or cutoff distances, which make the energy diverge.

My questions are as follows, and I would appreciate for your answers. Please also let me know if my questions need to be furthering specified.

[-in input file]

  1. unit: I certainly want to keep all the units reduced in the system, and it looks like the default lj fits most my needs. Is it correct?
  2. temperature: I need to assign temperature in pair_style and in velocity. As the example on the manual p.828, the temperature 1 is obviously in reduced unit, but in the example on p.1085 the velocity command is using Kelvin. Is it correct?
  3. pairwise cutoff: In the manual p.828-829, there are “cutoff” both in pair_style and in pair_coeff. I think the one in pari_coeff means the Rc of DPD (but in what unit?), then what is the cutoff in pair_style? Should it be the same as the one in pair_coeff, since we don’t have force without overlapping?
  4. neighbor & neibor modify: I manually tuned the values based on the suggestions from the error messages. I would like to know whether there is a specific way to decide the variables used in these two commands.
    [-in data file]
  5. For atom part, what I understand is “number, molecule number, atom type, coordinates x, y & z.” I learned from the micelle example in LAMMPS package that one can put all zero for the solvent molecules. I have 270 surfactant molecules, which I named it from 1 to 270 and each has 6 atoms. I also have another 79380 “solvent” molecules, which are named all in zero. Is it correct?

[input file]

This is LAMMPS input script specifies a DPD simulation of C8E8 (T2H4) in water

Parameters are assigned in in.c8e8, and configuration is in data.c8e8

File last edited 11/14/2012

dimension 3

neighbor 0.3 bin
neigh_modify one 10000 delay 3
communicate single vel yes

atom_style bond
pair_style dpd 1.0 30. 34387

read_data data.c8e8
bond_style harmonic

velocity all create 298.2 2349852

define masses and interaction coefficient

bead type:: 1: tail 2: head 3: water

mass * 1
pair_coeff * * 106.5 1.5 7.1
pair_coeff 1 3 126.1 1.5 7.1
pair_coeff 1 2 113.0 1.5 7.1
pair_coeff 2 3 107.5 1.5 7.1

define bond types and coefficients

bond type:: 1:c4-c4, 2:c4-eo2, 3:eo2-eo2, 4:c4-c4-eo2, 5:c4-eo2-eo2, 6:eo2-eo2-eo2

bond_coeff 1 80.0 0.8
bond_coeff 2 64.0 1.0
bond_coeff 3 53.3 1.2
bond_coeff 4 32.3 1.2
bond_coeff 5 29.1 2.2
bond_coeff 6 26.7 2.4

specify simulation parameters

timestep 0.01
thermo 1

first equilibrate the initial condition

fix 1 all nve
run 500

[data file]
LAMMPS 3d micelle c8e8 data file

81000 atoms
2430 bonds
0 angles
0 dihedrals
0 impropers

3 atom types
6 bond types
0 angle types
0 dihedral types
0 improper types

0.0000000E+00 30.00000 xlo xhi
0.0000000E+00 30.00000 ylo yhi
0.0000000E+00 30.00000 zlo zhi

Masses

1 1.000000
2 1.000000
3 1.000000

Atoms

1 1 1 9.433 8.694 5.254
2 1 1 9.487 8.082 5.808
3 1 2 9.360 8.882 6.348
4 1 2 8.423 9.582 6.641
5 1 2 7.637 10.584 6.149
6 1 2 6.451 10.579 5.762

81000 …

Bonds

Bonds

1 1 1 2
2 4 1 3
3 2 2 3
4 5 2 4
5 3 3 4
6 6 3 5


2430 …

[errors]

LAMMPS (28 Aug 2012)
Scanning data file …
2 = max bonds/atom
Reading data file …
orthogonal box = (0 0 0) to (30 30 30)
2 by 2 by 2 MPI processor grid
81000 atoms
2430 bonds
Finding 1-2 1-3 1-4 neighbors …
4 = max # of 1-2 neighbors
8 = max # of 1-3 neighbors
23 = max # of 1-4 neighbors
5 = max # of special neighbors
Setting up run …
Memory usage per processor = 129.433 Mbytes
Step Temp E_pair E_mol TotEng Press
0 298.2 85036.554 0.032202294 85483.881 256595.23
1 1.8911299 85033.665 0.033103112 85036.535 255708.16
2 5.9758557 85038.894 0.033937598 85047.891 255743.6
3 41.145388 85086.615 0.039222409 85148.372 256147.47
4 414.38357 85592.155 0.076058462 86213.799 259721.24
5 3877.0471 90182.595 0.34279814 95998.436 289944.43
6 12538.684 99473.059 0.95076323 118281.8 300376.21
7 65645.291 47130.775 11.092103 145608.59 190320.21
8 46377.024 30448.454 51.14401 100064.27 188363.08
ERROR on proc 6: Bond atoms 369 371 missing on proc 6 at step 9 (neigh_bond.cpp:55)
ERROR on proc 7: Bond atoms 543 544 missing on proc 7 at step 9 (neigh_bond.cpp:55)
ERROR on proc 0: Bond atoms 904 906 missing on proc 0 at step 9 (neigh_bond.cpp:55)
ERROR on proc 3: Bond atoms 574 576 missing on proc 3 at step 9 (neigh_bond.cpp:55)
4 additional processes aborted (not shown)

Several of your Qs are about units. The rule for units
is simple. Whichever units you use (via the units command),
all inputs and outputs to/from a simulation are in those units.

More comments below.

Steve

I am initiating a studying case of micellization of a nonionic surfactant in
water. The simulations were terminated in a few steps due to the divergent
energy, and I figured it could be my misunderstandings for some parameters
introduced in the LAMMPS manual. What I did is addressed as follows. Input
file and the data file are described in the bottom.

The system contains two molecule types, surfactant and water. A surfactant
is modeled by 6 beads (1-1-2-2-2-2) containing two bead types 1 and 2, and
water is modeled as bead type 3. Bead size (Rc) is unified as 7.1 nm, and
all length units in the system are reduced by this value (box length, bead
coordinates, bond length, etc.) Bead mass is unified as 1. The simulation
box is 30x30x30 (Rc^3). Surfactant bond types (two body for connection and
three body for rigidity) are assigned using harmonic bonds. Number of beads
is 81000, projecting to the density equals to 3 in the box. Gamma
(force/velocity unit) in DPD is set to 1.5. Initial coordinates are taken
from the final configuration of my previous calculations under same
conditions. After performing parallel simulations using 8 cores, I had error
message "Bond atoms %d %d missing on proc %d at step %ld." I think that I
must set the wrong pairwise parameters or cutoff distances, which make the
energy diverge.

That error is usually due to bad dynamics. See Section_error in the manual.
If a bond blows apart, then
one processor can't find the atom it needs b/c it is too far away on another
proc. Run on one proc until you are sure the dynamics is good (viz, thermo
output). You could also get that error if your DPD cutoff is shorter than the
max bond length. In that case, see the communicate cutoff command to
grab atoms from further away than the pair cutoff + neighbor skin. See
the communicate doc page for details.

My questions are as follows, and I would appreciate for your answers. Please
also let me know if my questions need to be furthering specified.

[-in input file]
1. unit: I certainly want to keep all the units reduced in the system, and
it looks like the default lj fits most my needs. Is it correct?

that's up to you

2. temperature: I need to assign temperature in pair_style and in velocity.
As the example on the manual p.828, the temperature 1 is obviously in
reduced unit, but in the example on p.1085 the velocity command is using
Kelvin. Is it correct?

the T is in whatever units you select, could be reduced, could be Kelvin

3. pairwise cutoff: In the manual p.828-829, there are "cutoff" both in
pair_style and in pair_coeff. I think the one in pari_coeff means the Rc of
DPD (but in what unit?), then what is the cutoff in pair_style? Should it
be the same as the one in pair_coeff, since we don't have force without
overlapping?

The DPD doc page explains this. The style cutoff is the general default
cutoff for all atom type pairs. The pair coeff cutoff can override it
for a specific pair. If you don't want to override it, you don't need to
specify it.

Again, the doc page says the sigma is in "distance" units. What that
means depends on what units you define. If it's LJ, then it's sigma.
If it's real units, it's Angstroms, etc.

4. neighbor & neibor modify: I manually tuned the values based on the
suggestions from the error messages. I would like to know whether there is a
specific way to decide the variables used in these two commands.
[-in data file]

They are settings. The defaults are usually fine. If not then set
them to something more appropriate for your system. Ask a more
specific question if the doc page doesn't explain them sufficiently.

5. For atom part, what I understand is "number, molecule number, atom type,
coordinates x, y & z." I learned from the micelle example in LAMMPS package
that one can put all zero for the solvent molecules. I have 270 surfactant
molecules, which I named it from 1 to 270 and each has 6 atoms. I also have
another 79380 "solvent" molecules, which are named all in zero. Is it
correct?

You can assign atoms any molecule ID you want (>= 0 I think). So there
is no "correct" value. What you describe is fine.