[lammps-users] instability/bug in DPD-like potential

Hi all.

I have written a pair potential which does local damping (Kelvin damping a la the fracture simulations by Holian et. al.). It is implemented in an analogous way to pair_dpd.cpp, but damps the relative velocity vector rather than just the radial component. It seems to work fine for the most part, but there is some bizarre behavior which I’m hoping someone has some thoughts on. It works well enough otherwise that I’ve been using it for a couple of months without noticing this bug.

In order to see the strange behavior, it seems like I need all of the following 3 things:

  1. a multi-processor simulation
  2. a volume/rescale fix on one direction and a zero-pressure nph on the other (simulations are 2D) (volume/rescale in both directions won’t produce the bug… maybe non-zero confining pressures will produce the bug too, but I haven’t tried)
  3. my home-brewed Pair class which implements the Kelvin dissipation

The behavior is as follows: eventually, one of the pairs of particles which lies near a CPU boundary gets locked into a very tight orbit (few tenths of a particle diameter) at large velocities. The interaction potential is lj_smooth, and max velocities are order 10 in LJ units… period is only a few timesteps, and timesteps are 0.0025 in LJ units. The displacements seem to stay along a single ray with some slight drift, and the velocity just keeps changing sign but not direction. The particles never really get that far apart, and one would only notice something pathological going on by looking at, e.g., kinetic energy, since I don’t add the viscous contribution in to the virial.

I call it a “bug” rather than a numerical instability, since I can only get it to happen when I do a multi-processor run.

Here’s a relevant snippet from my pair class which comes after computing the usual pairwise forces:

-------------------- pair_lj_smooth_kelvin.cpp -------------
… snip …
// Now put on the drag forces:
double delvx = v[i][0] - v[j][0];
double delvy = v[i][1] - v[j][1];
double delvz = v[i][2] - v[j][2];

f[i][0] -= delvxvisc[itype][jtype];
f[i][1] -= delvy
f[i][2] -= delvz*visc[itype][jtype];
… snip …

Is the time evolution of the 2 particles in an orbit correct?
I.e. treated by themselves are their forces and velocities and positions
being timestepped correctly? If not, then print out what's contributing
to their motion and see what's going wrong.

I assume you are using atom_style dpd which updates ghost atom
velocities each step.


No. I wasn’t. This makes perfect sense and seems to fix the problem. Thanks!

Actually, check that…

problems still creep up (at later times) with atom_style = hybrid molecular dpd. Presumably it’s the same kind of problem. I’ll dump even more detailed debugging info about the bad guys and see what’s going on with the forces and the timestepping…

I also changed the Kelvin dissipation, so that there is a smooth cutoff like in pair_dpd. Doesn’t seem to help.