Handling weak bond potentials

LAMMPS LAMMPS Beginners
1

I am trying to reproduce the simulations from the paper cited in all the LAMMPS DPD documentations: Groot and Warren, J Chem Phys, 107: 4423-4435 (1997). DOI: 10.1063/1.474784.

I am now having troubles simulating the bonded systems (simple polymer melt of a few chains) they describe.
In principle, it is just a simple

bond_style harmonic 
bond_coeff * 1.0 0.0

, or alternatively, a

bond_style table linear 20
bond_coeff * linear_bond_coeff_2.table LIN

and some bonds in the initial structure (atom_style hybrid bond edpd), together with the pair style:

special_bonds lj 1.0 1.0 1.0 
special_bonds coul 1.0 1.0 1.0
pair_style dpd 1.0 1.0 ${seed}
pair_coeff * * 25 4.5

The issue is now, that this bond force is very weak. While LAMMPS happily runs this simulation, and I can reproduce e.g. the mean bond length, I am getting problems when I want to investigate results of properties that require the absolute coordinates, such as the radius of gyration or the mean square displacement, as, because of the weak bond force, some bonds reach a length > 0.5 * the box length, in which case the particle is pulled to the next periodic image. A system more than 100 times the size reported by Groot Warren has still sometimes shows this behaviour.

One may argue, that the bond forces should be computed in absolute terms. As I understand it, LAMMPS doesn’t, and behaves the way it does (computing the bond forces in a PBC corrected manner), in order to easily support infinite lattices and networks.

In literature, one can find enough references of people using infinite networks with this weak-bonded DPD potential, leading me to belive that there are “simple” ways to support both.
My question: are there “simple” ways also to introduce support of such weak bond terms into LAMMPS? Or do I have to find a different software package if I want to use such a weak bond potential with such strong pair forces/implement my own?

2

Not simple at all. Your experiments only work with a domain decomposition code when running on a single processor. Once you run in parallel, you would sooner or later get a “bond atom missing” error. I would seriously doubt that there is any general purpose MD code that handles the situation your are describing the way you want it to be handled. Doing this in parallel is next to impossible with domain decomposition. LAMMPS has a special purpose pair style list to work around the domain decomposition limitations (hence disguising the bond interaction as a pair style), but that implementation is still subject to minimum image conventions, so you cannot have a bond that stretches longer than half the box in any direction.

1 Like
3

Thank you for your answer. This is what I suspected.

However, I could not actually believe that the bonds would get as long as they did; I did not consider it realistic and would not understand, why literature would use this force field like that.
Upon further review, I found that the default neighbour list settings were once again hindering; adjusting those leads to maximum bond lengths that are still quite long, but at least do not reach > L/2 anymore for sensibly sized systems.