Failure to Replicate DPD Simulation

I am trying to replicate the DPD simulation in LAMMPS. With that, I mean the simulation leading up to Figure 1 in the paper cited in all the LAMMPS DPD documentations: Groot and Warren, J Chem Phys, 107: 4423-4435 (1997). DOI: 10.1063/1.474784.

Now, I can actually replicate the results they get, but only shape and tendency wise, the data I find is offset on the x-axis compared to the paper’s, and it is not clear to me, why this would be the case.

This is the plot I try to replicate:

and this is how it currently looks like:

As you can see, the x-axis is very different.
My LAMMPS input looks like this (standalone, you should be able to run like this):

# LAMMPS input file: fig. 1 replica of Groot, Warren, 1997
units           lj
atom_style tdpd 1
variable kb equal 1.0
comm_modify vel yes

variable seed equal 2570894

# create the box and atoms with rho = 4
region thebox block 0 10 0 10 0 10
create_box 2 thebox 
create_atoms 1 random 4000 ${seed} thebox
mass * 1

variable sigma equal 3
variable T equal 1. 
variable gamma equal 0.5*${sigma}*${sigma}/(${kb}*${T})
variable A equal 25
variable cutoff equal 1.0

pair_style dpd ${T} ${cutoff} ${seed}
pair_coeff * * ${A} ${gamma} ${cutoff}

fix 1 all mvv/tdpd 1.0
# run_style verlet 

run 100

# output we are interested in
compute 3 all temp/com
thermo_style    custom step temp dt press c_3 vol
thermo 1

# dump 1 all custom 1000 ./output/fig-1.dump.txt id fx fy fz

# actual simulation
variable timeStp index 0.0015 0.0025 0.0035 0.005 0.007 0.01 0.015 0.02 0.0025 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.2 0.25 # it will loose atoms at some point, but we don't care
label loop0
# adjust the time-step
timestep ${timeStp}

# re-run the simulation with the new time-step
run 20000

next timeStp
jump SELF loop0

and my LAMMPS version is LAMMPS (23 Jun 2022), as installed on MacOS using Homebrew.

I would expect to be able to use LAMMPS in accordance with the documentation to reproduce said paper, but it is possible that I am misunderstanding something.
Some observations:

  • I require much more than the 200 steps Groot & Warren report having done to get a stable result. This would hint, that either they did not have a stable result, or the random forces I get are too strong.
  • I would consider the general behaviour to be reproduced, ignoring the offset on the x-axis.

I would kindly ask for some feedback on what could be reasons for this replication to fail, or if anyone notices any issues in the code, I would like to understand why my simulations break down at a different time-step.

1 Like

Why atom style tdpd and fix mvv/tdpd?

The docs for the fix say:

Fix mvv/tdpd updates the per-atom chemical concentration, in addition to position and velocity, and must be used with the pair_style tdpd command.

I would expect that you use the default atomic pair style and fix nve.

Thanks for your answer, dear @akohlmey.

I want to use one of the mvv/dpd fixes for integration in order to actually reproduce the different lambdas, as this is not possible with fix nve.
Now, to use a mvv/dpd type of integrator fix, I need to have an atom_style with an attribute vest. Not even atom_stlye dpd does that, so I simply switched to tdpd, which is one of the few that does, as I don’t mind a constant concentration.
With this, I can also use mvv/dpd instead of mvv/tdpd, the results are the same, but mvv/tdpd seemed to me better in the code I post due to consistency.

1 Like

Atom style dpd belongs to a very different package and thus has no relation to those fixes in the DPD-MESO package. Any of the atom styles in the package should work. edpd seems to be the minimal choice there.

At any rate since you do see differences, the obvious way to proceed is to eliminate unknowns and using the DPD-MESO fixes is far less tested that fix nve. So even if you are locked into simulating a subset or options, the first thing to do would be to check if the result with fix nve is consistent with using lambda = 0.5 and one of the integration fixes in DPD-MESO.

Unfortunately, I cannot provide any more specific advice, since DPD is outside my area of (scientific) expertise.

Thank you for your answers. I will keep updating this post nonetheless with my progress for future searchers.

Using fix nve is indeed consistent with the integration fixes from DPD-MESO, at least where ca. dt >= 0.02.

There are only two components to this simulation: an integrator, and a force-field, which is only composed of a pair style. Due to the consistency with the nve integrator, I would expect that a search for the error in the pair style would be more sensible now, yet changing the two parameters I have (gamma, T) away from the reported values leads to shapes that are no longer consistent with the results reported in 1997.

1 Like

I have just one more remark. Pair style dpd has been around in LAMMPS for a very long time. While we had some bugs in the code that went unnoticed for a long time, those were usually for corner cases or rarely used features. DPD does not fall info those categories, so it is more likely there are some differences in your input or system setup compared to what the simulations in the reference were using. Perhaps you misread something or there is a unit conversion or convention issue. This is all difficult to assess from remote and not really a LAMMPS issue as well.

After implementing the whole algorithm myself, separately, to get a third opinion on this whole thing, I was able to get some hints that the results by Groot/Warren are indeed correct.

I went through the relevant code in both codebases, LAMMPS and mine, but did not find notable differences. Finally, I recompiled LAMMPS, with which I was suddenly able to reproduce the results from Groot/Warren, with the same input file. As I have both versions of LAMMPS, still, I can really run the same input to get two different outputs. With the NVE as well as the mvv/dpd.

I could not find a relevant change in the LAMMPS source code yet. No idea if it is related to a parameter not being read correctly, or some strange reason like that, but anyway, while it would help for my peace of mind to know what the issue was back then, at least I can progress now.

1 Like

After spending another many hours trying to figure out what the difference between the two LAMMPS versions is, I resorted to recompiling LAMMPS for different commits to enclose the moment when the change happened. I tried with a binary search, but due to many merges from other people, it seems that my script was not able to get the history to line up linearly and there are working and not-working versions alternating.

Now, after having compiled 75 different commits, I finally found the origin: in 082defa862fc2db66e15ce66713359cb868f878e, the default of the neighbour list update was changed. With the old default, the temperature with dt = 0.03 was 10 times higher than with the new default.
This is reproducible with neigh_modify.

It took a while, but I am happy that I figured it out. This should restore my trust in LAMMPS (or my LAMMPS skills, respectively).

1 Like

Is there a way, dear @akohlmey , how I can mark this topic as resolved or something like that?

@GenieTim I am happy to see that you figured out what was causing the trouble.

As for marking the topic as solved, I have not noticed an option to do this here. I’ve seen a “solved” marker in some other categories, but have no clue how this was configured or who needs to set this up.

What can certainly be done is to “close” the topic (but I am not certain if that is a good idea, since somebody else may have a similar issue and a followup question). So your statement of having resolved the issue and the description of how is probably the best that can be done.

1 Like