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
(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)?
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.