Questions about abnormal temperature fluctuation

Hello. Recently I encounter some problem in running LAMMPS. The system is modeling micelle of a linear surfactant using DPD. I ran it on LAMMPS (not parallel) comparing to the results from group code, using the same parameters. An additional 1-3 bond is applied as surfactant rigidity. While running in LAMMPS w this 1-3 bond, there is a sign that system is heating up. In nve, the temperature will keep climbing even to several thousands in reduced unit. If it is in nvt, the reduced temperature will fluctuating in a range around 1 to 3. Such level of fluctuation was never observed while using group code.

I have checked my inputs very carefully and over and over with LAMMPS manual, but no error is found by myself. I have also checked my input scripts with another experienced LAMMPS users, but the inputs looks also fine to her. I was wondering whether I could submit my data and input files to the forum, and my questions and systems will be furthering described. It is a 3000-atom testing system, so the files won’t be huge.

If it is inappropriate to do so, would someone please provide a hint that where the problem could be and what kind of test should I do. Suggestions are appreciated.

sure - you can submit them - but it's not likely
someone will debug it carefully for you. If you have
another code to compare to, I suggest you compare
snapshots and see if forces and energy are identical.

Steve

Thank you Steve for the suggestions. I am comparing two codes and the output trajectories. I don’t expect someone will resolve the problem for me. I just want to confirm that the situation I have doesn’t come from some simple mistakes.

Problem descriptions:
Using the dissipative particle dynamic simulation, I want to simulate the motion of surfactant (type 1 for hydrophobic tail and type 2 for hydrophilic head) molecules in the water solvent (type 3 atom). Box size is 20*Rc per side, and there are only 80 6-bead surfactant in 23520 water bead (for reduce density=3). Besides the spring forces for surfactant bead connections, I employed another type of spring force on 1-3 body, as the rigidity of the molecules. However, once the 1-3 body harmonic bonds are added, there is a big temperature fluctuation in the simulations. Please check “out.c8e8-2” for thermo w/o 1-3 bonds, and out.c8e8-3 for thermo w 1-3 bonds.

One can reproduce the case by using attached input and data files, and 1-3 bonds can be added by changing the comments line in the input file. The default simulation step is short in the input file, but I have checked that the fluctuation won’t go away with long run. Please kindly let me know whether the situation is normal, or it is some error resulted from improper inputs. Please let me know the problem needs to be furthering described. Thank you.

data.c8e8 (715 KB)

in.c8e8 (1.11 KB)

out.c8e8-2 (8.61 KB)

out.c8e8-3 (8.05 KB)

A couple of comments.

a) I would run NVE to see what the dynamics do, before using
a thermostat.

b) The pair coeff value of 100+ seems large for the 1st param
I'm used to seeing DPD models with a value around 1, like
everything else in LJ units.

c) Have you verified that the 2 codes give the same energies/forces
(pressures) for a few snapshots of your system. If so, I find
it hard to believe that the dynamics would be very different,
especially for NVE.

Steve

Hi Steve,

I am working on ©, but I really want to know your further comments on (a) and (b). Thank you.

(a) I actually use NVE in the beginning. Temperature using nve with the additional 1-3 bonds will be increasing during simulation, and that’s one reason I switch to thermostat. In the case without 1-3 bonds, temperature and other thermo output in nve are similar to that in nvt, and that’s what I expect for DPD. It looks like some heat is generating with the additional 1-3 bonds, and I just don’t know where it comes from.

(b) I would appreciate if you can explain the “force unit” in DPD in more details. I am worry about I may misunderstood according to what You said that the common value for repulsion parameter is one in LJ unit. In the beginning, I used Dr. Malibaum’s post as a input template (http://lammps.sandia.gov/threads/msg05160.html), the pair coefficient was the “most common” value 25. Using the parameters 25, I was able to reproduce the pressure in Groot and Warren’s paper. The repulsion parameters are larger in some papers, such as 78~104 in Groot and Rabone’s work (Biophysical Journal Volume 81 August 2001 725–736). In LAMMPS unit page, force needs to be rescaled by epsilon/sigma. In pair_dpd.cpp, I didn’t see any rescaling of “a0”. However, do I somehow need to rescale the repulsion parameters from those values on the papers (say 25 on G&W and ~100 on G&R)?

Comments below.

Steve

Hi Steve,

I am working on (c), but I really want to know your further comments on (a)
and (b). Thank you.

(a) I actually use NVE in the beginning. Temperature using nve with the
additional 1-3 bonds will be increasing during simulation, and that's one
reason I switch to thermostat. In the case without 1-3 bonds, temperature
and other thermo output in nve are similar to that in nvt, and that's what I
expect for DPD. It looks like some heat is generating with the additional
1-3 bonds, and I just don't know where it comes from.

It is puzzling to me also. Since you say you have another code to
compare to, I am saying you should compare the 2 codes, with and without
the added 1-3 terms. Assuming you have chosen parameters (e.g.
the bond strength and bond length) correctly for the 1-3 terms, the only
thing I can imagine would cause a rise in temperature is if your initial
config is far from equilibrium (with respect to the 1-3 bond length). For
example if you want the chains to be nearly linear (e.g. the 1-3 bond
length R0 is roughly 2x the 1-2 bond length R0), but you start with
chains that are all coiled up, then the temperature could rise, b/c there
is a lot of initial energy in those 1-3 bonds.

(b) I would appreciate if you can explain the "force unit" in DPD in more
details. I am worry about I may misunderstood according to what You said
that the common value for repulsion parameter is one in LJ unit. In the
beginning, I used Dr. Malibaum's post as a input template
(LAMMPS Molecular Dynamics Simulator), the pair coefficient was
the "most common" value 25. Using the parameters 25, I was able to reproduce
the pressure in Groot and Warren's paper. The repulsion parameters are
larger in some papers, such as 78~104 in Groot and Rabone's work
(Biophysical Journal Volume 81 August 2001 725–736). In LAMMPS unit page,
force needs to be rescaled by epsilon/sigma. In pair_dpd.cpp, I didn't see
any rescaling of "a0". However, do I somehow need to rescale the repulsion
parameters from those values on the papers (say 25 on G&W and ~100 on G&R)?

Don't think of units being specific to DPD. In LAMMPS, there
is just LJ units (or real or metal etc). All the DPD doc page
is saying is that a0 is used as a pre-factor on a force term,
so it needs to be in force units (LJ or metal or real, etc). All
I was saying is that in LJ units it is common for every parameter
to be close to 1.0 (sigma, epsilon, temperature, energy, etc).
A particular param might be 0.5 or 2, but 100 would be pretty
unusual. Hence I'm suspicious that you are running the model
you expect. Again, comparing to another code that does what
you want would be a good way to check this.

Hi Steve,

Thank you for the reply. I checked the snapshots from two codes, and they yielded different results w & w/o 1-3 bonds. I am not comparing the codes (yet), but only exploring the case in LAMMPS. You are probably right about that I am not running the things I wanted in LAMMPS, although I still don’t know where the problem could be. I will try another configuration with bond length close to equilibrium length, and see whether there is still a fluctuation. I understand that there will be a heat-up if I start with length longer than equil. length, but the temperature should goes down with the bond length approaches to its equil. length with simulation going. Isn’t it?

Also, I am confused by what you said about force prefactor A. I understand the concept of lj unit is to keep everything around the order of unity. Are you suggesting me to rescale my force parameters by dividing some trivial epsilon number to make every parameter close to 1? For example, if I divide all my around-100 A by 100, I can have all repulsion parameters near 1. Nevertheless, gamma (“force”/velocity units) in pair_style dpd also needs to be divided by 100, and it will generates a very small number. As I said, I am able to reproduce the pressure in Groot and Warren’s 97’ paper by using A=25 (which is also >>1) in lj unit, and that’s exactly why I think the A in LAMMPS is what it is in most DPD papers. A~100 is actually pretty common since Groot and Rabone’s work. May I have more your thoughts on this? Thank you.

Sincerely,
Ming-Tsung

I'm not very familiar with DPD models, so maybe 25 or 100 is fine.
It's just large for LJ units.

Steve