Getting the DPD thermostat right

Dear all,

I have spent a considerable amount of time trying to make DPD with the thermostat work correctly, yet my results are clearly wrong. Can you please help?

Reading the manual, I understand that there are in general two options for thermostatting, starting from temperature T, the cutoff, interaction parameter a, friction gamma, and the seed:

pair_style dpd T cutoff seed
pair_coeff * * a gamma
fix 1 all nve
fix 2 all langevin T T gamma seed

pair_style hybrid/overlay dpd T cutoff seed dpd/tstat T T cutoff seed
pair_coeff * * dpd a gamma
pair_coeff * * dpd/tstat gamma
fix 1 all nve

As far as I understand, the first option makes the particles interact with an implicit solvent and not each other and so is not a proper DPD. Hence, only the second option remains.

Trying the second option on a pure DPD liquid, the temperature goes up to 3.5 and does not settle (it should settle at 1.0), and the RDF is totally incorrect close to r = 0.

Can I please ask what exactly wrong with my input, and how to get it right? I tried to follow the manual as closely as it gets. Attached I am sending the data and input file together with the RDF.

Thank you in advance.

Kind regards,

rdf_1.png

binmixt.in (1.37 KB)

binmixt.data (98.6 KB)

Dear all,

I have spent a considerable amount of time trying to make DPD with the
thermostat work correctly, yet my results are clearly wrong. Can you please
help?

Reading the manual <http://lammps.sandia.gov/doc/pair_dpd.html>, I
understand that there are in general two options for thermostatting,
starting from temperature T, the cutoff, interaction parameter a, friction
gamma, and the seed:

1.
pair_style dpd T cutoff seed
pair_coeff * * a gamma
fix 1 all nve
fix 2 all langevin T T gamma seed

2.
pair_style hybrid/overlay dpd T cutoff seed dpd/tstat T T cutoff seed
pair_coeff * * dpd a gamma
pair_coeff * * dpd/tstat gamma
fix 1 all nve

​both of these inputs make no sense.
the effect of pair style dpd/tstat is *already included* in pair style
dpd.​ so there is no point in using hybrid/overlay. the reason for
dpd/tstat is to apply it to other potentials.

similarly, it is not correct to use fix langevin in combination with pair
style dpd. you are applying a thermostat twice with conflicting approaches.
that is likely to be bogus.

​axel.​

​​

Dear Axel,

thanks a lot for the clarification. So I understand that the input below fully describes the DPD in LAMMPS:

pair_style dpd T rc seed
pair_coeff * * a gamma rc
fix 1 all nve

I am not sure if I am the only one, but I think if would be very useful if these lines were stated in the DPD section of the manual.

To improve my understanding I am now trying to reproduce a more complex system – phase separation of a symmetric diblock copolymer, as described in this paper. Polymerisation N = 10, like bead interaction 25 and cross interaction 50. In theory I should be easily getting a lamellar phase after as few as 20k steps, but in practice I am see only a disordered phase. However, I was able to get this phase with the DL_MESO package. I am sending the screenshots and input file attached.

Can you please suggest what I am doing wrong? I have tried everything I could think of.

Thank you in advance.

Kind regards,

diblock.in (1.26 KB)

vmdscene_dlms.png

vmdscene_lmp.png

Dear Axel,

thanks a lot for the clarification. So I understand that the input below
fully describes the DPD in LAMMPS:

pair_style dpd T rc seed
pair_coeff * * a gamma rc
fix 1 all nve

I am not sure if I am the only one, but I think if would be very useful if
these lines were stated in the DPD section of the manual.

​you are the first person in almost 10 years that i recall being confused
in this way by what is in the docs about the DPD pair style. i think the
description is very clear that dpd/tstat is basically dpd​ with some stuff
left out to make it a thermostat only. so that makes your second option
illogical. and the first option you used is strange as well because you
already include a temperature parameter in the pair style command. add to
this the pointer to a publication explaining the entire model the style is
implementing, i don't see how this is a significant improvement. not to
mention, that fix nve is not the only integrator compatible with this pair
style, fix nph should also work.

To improve my understanding I am now trying to reproduce a more complex
system -- phase separation of a symmetric diblock copolymer, as described in
this paper <http://aip.scitation.org/doi/abs/10.1063/1.476300>.
Polymerisation N = 10, like bead interaction 25 and cross interaction 50.
In theory I should be easily getting a lamellar phase after as few as 20k
steps, but in practice I am see only a disordered phase. However, I was
able to get this phase with the DL_MESO package. I am sending the
screenshots and input file attached.

Can you please suggest what I am doing wrong? I have tried everything I
could think of.

​sorry, i don't have the time to review people's work and correct it. if at
all, this would be the task of your adviser/supervisor.

axel.

Dear Axel,

I am writing after a year, having spent long hours working with LAMMPS. I have two points.

  1. You said in your previous email:
    “​you are the first person in almost 10 years that i recall being confused in this way by what is in the docs about the DPD pair style.
    Given you are well-known for your passive-aggressive responses, I can only guess how keen people are to write if they are confused about something in LAMMPS. There might be some bias, like there were no complaints about human right abuses in Stalinist Russia.
    More importantly, in December 2015 I wrote an email about my confusion with thermostatting DPD, and you responded: “i have to admit that the pair style dpd documentation is lacking a bit in this regard.” I find these two claims inconsistent. Why is it surprising that I am confused, if you yourself admit that the documentation is lacking?

Hence, I suggest adding the attached file dpd.in to the DPD page of the manual as a standard example of how to simulate water, or any single-component liquid. At least the critical lines of it. (The use of SI units was advised to me by Sergey Litvinov, I think LJ units would be as good). I am sure many readers would appreciate that.

  1. A year ago, I mentioned in this thread that simulating polymers in LAMMPS fais to produce correct results, according to standard literature (Groot, JCP, 1998). Recently, I have spent a lot of time investigating this problem, which is the reason I am writing again. I have thoroughly compared the results from LAMMPS with DL_MESO (denoting DLMS), another package with implemented DPD.

I have done the following runs, with like species interaction parameters “a” and unlike “a+da”, timestep dt=0.05, 10k steps:

  • DLMS, binary mixture, a=25, da=10, correct result, phase separation visible in VMD,
  • LMP, binary mixture, a=25, da=10, correct result, phase separation (attached input file binmixt.in),
  • DLMS, diblock copolymer, a=25, da=10, correct result, lamellar phase established,
  • LMP, diblock copolymer, a=25, da=10, wrong result, no phase separation (attached input file diblock.in).

Each time I monitored the radial distribution function (RDF), which can reveal many problems.

It is visible that the RDF from LAMMPS (rdf_lmp_diblock.png) produces a massive peak at r=0, which is a very strange behaviour. Compare with DLMS (rdf_dlms_diblock.png).

I tested a limiting case, setting “bond_style harmonic” to 0 0. The polymer melt with freely moving beads is now essentially a binary mixture (for correct behaviour, see rdf_lmp_binmixt.png). In this case VMD showed a clear phase separation, but the RDF is again wrong, with a peak at r=0 (see rdf_lmp_diblock_zero_spring.png).

At this point, the only difference in the input files between binary mixture (binmixt.in) and diblock copolymer with zero spring constant is “atom_style atomic” vs “atom_style molecular”/“atom_style bond”. Clearly, redefining the atom_style changes physics. Here I am really confused and not sure how I can proceed further to squeeze the correct behaviour from LAMMPS, which is the reason I am writing again.

I believe this is strong enough case for you to take notice. I fully realise that LAMMPS is a free open-source code with no warranty and you are in no way obliged to solve my problems. But I think you might be interested.

Thanks in advance for consideration.

Peter

dpd.in (972 Bytes)

binmixt.in (1.16 KB)

diblock.in (1.24 KB)

rdf_dlms_diblock.png

rdf_lmp_binmixt.png

rdf_lmp_diblock.png

rdf_lmp_diblock_zero_spring.png

Dear Axel,

I am writing after a year, having spent long hours working with LAMMPS. I have two points.

  1. You said in your previous email:
    “​you are the first person in almost 10 years that i recall being confused in this way by what is in the docs about the DPD pair style.
    Given you are well-known for your passive-aggressive responses, I can only guess how keen people are to write if they are confused about something in LAMMPS. There might be some bias, like there were no complaints about human right abuses in Stalinist Russia.
    More importantly, in December 2015 I wrote an email about my confusion with thermostatting DPD, and you responded: “i have to admit that the pair style dpd documentation is lacking a bit in this regard.” I find these two claims inconsistent. Why is it surprising that I am confused, if you yourself admit that the documentation is lacking?

i don’t think that those two statements are in contradiction. lots of people have used LAMMPS DPD (and there are now lots of more elaborate versions of DPD available as well). since nobody there has not been any complaint for such a long time, i believe i can assume, that for all of those people the docs were sufficient.

​also, in those two e-mails your are confused about two very different things. and the behavior you are describing in your december 2015 e-mail is not unusual. the questions you ask are those of a person with limited experience in MD fundamentals. ​

Hence, I suggest adding the attached file dpd.in to the DPD page of the manual as a standard example of how to simulate water, or any single-component liquid. At least the critical lines of it. (The use of SI units was advised to me by Sergey Litvinov, I think LJ units would be as good). I am sure many readers would appreciate that.

  1. A year ago, I mentioned in this thread that simulating polymers in LAMMPS fais to produce correct results, according to standard literature (Groot, JCP, 1998). Recently, I have spent a lot of time investigating this problem, which is the reason I am writing again. I have thoroughly compared the results from LAMMPS with DL_MESO (denoting DLMS), another package with implemented DPD.

I have done the following runs, with like species interaction parameters “a” and unlike “a+da”, timestep dt=0.05, 10k steps:

  • DLMS, binary mixture, a=25, da=10, correct result, phase separation visible in VMD,
  • LMP, binary mixture, a=25, da=10, correct result, phase separation (attached input file binmixt.in),
  • DLMS, diblock copolymer, a=25, da=10, correct result, lamellar phase established,
  • LMP, diblock copolymer, a=25, da=10, wrong result, no phase separation (attached input file diblock.in).

Each time I monitored the radial distribution function (RDF), which can reveal many problems.

It is visible that the RDF from LAMMPS (rdf_lmp_diblock.png) produces a massive peak at r=0, which is a very strange behaviour. Compare with DLMS (rdf_dlms_diblock.png).

I tested a limiting case, setting “bond_style harmonic” to 0 0. The polymer melt with freely moving beads is now essentially a binary mixture (for correct behaviour, see rdf_lmp_binmixt.png). In this case VMD showed a clear phase separation, but the RDF is again wrong, with a peak at r=0 (see rdf_lmp_diblock_zero_spring.png).

At this point, the only difference in the input files between binary mixture (binmixt.in) and diblock copolymer with zero spring constant is “atom_style atomic” vs “atom_style molecular”/“atom_style bond”. Clearly, redefining the atom_style changes physics. Here I am really confused and not sure how I can proceed further to squeeze the correct behaviour from LAMMPS, which is the reason I am writing again.

​having bonds with no force constant and having no bonds is resulting in different behavior due to “exclusions”. even if there are no forces from the bonds, the corresponding pairs of atoms are excluded from the neighbor list. you also have to adjust the exclusions settings via the “special_bonds” command as well to make those two settings to represent the exact same physics. again, you are showing lack of experience in ​MD fundamentals, and thus are not adding much credibility to your findings.

I believe this is strong enough case for you to take notice. I fully realise that LAMMPS is a free open-source code with no warranty and you are in no way obliged to solve my problems. But I think you might be interested.

​no, i am not.​ you have to do a better job to provide convincing evidence. the statistics are stacked against you. very experienced people have used DPD in LAMMPS and while it is not impossible, that there is a bug, it is extremely unlikely, and you are not providing any convincing evidence. if you get different results with a different software, it may just be, that that software has different default settings in parameters, that you don’t set explicitly (e.g. exclusions) and don’t check for.

​axel.​

this has nothing to do with the science, so i’ll put it in a separate response.

i don’t think passive-aggressive is a good description. “active aggressive” would be more fitting. or even better is “demanding” or “challenging” or “direct”

i am responding to e-mails on mailing lists for a very long time (not just on LAMMPS) and i have seen all kinds of people and all kinds and styles of questions. by responding the way i do, i aim to reduce wasting time (mostly mine). yes, it is more frustrating for somebody like you. no question. but also, there are many people asking questions and very few people answers, so that puts people answering at a disadvantage, especially when people come with the “i am confused, please tell me what is wrong with what i am doing”. it is not the job of this mailing list to teach you your craft. perhaps you can entice somebody else, that actually has practical experience with DPD (i don’t), but if you want to convince me, to take a second look you have to be more explicit that “i get difference results with a different code”. i don’t have the time to redo your work and study and compare in great detail both code’s documentation, build (small and trivial) test inputs and compare them with what i can compute manually and so on. this is all part of learning the craft of MD and that is somehting between you and your adviser. not you and me.

i know this doesn’t really solve your problem, but this is how i am and how i see the situation. you are free to ignore me and i will not press this issue any further. yet to make a case, that would make me become interested, you have to do a much better job and provide much more (technical detail) and show more convincing arguments. the fact, that you have spent a lot of time on this, is irrelevant. some people work on problems for a long time, usually because they could not find a way to understand them. only when you really understand a problem, it is easier to solve. the best strategy to solve complete/frustrating problems is to make them smaller and solve only parts. this is common sense and part of how research is conducted. however, it is not a topic that is suitable for a mailing list to discuss (technical) issues with using, programming or compiling/installing LAMMPS.

axel.