Dear all,
I am recently try to apply the DPD thermostat (pair dpd/tstat) on a simple Couette flow to control the temperature,
that’s where I have questions. For an equilibrium system, dpd controls the temperature via the dissipative and random]
force terms. When the system is hot, the relative velocities between particles is large and the dissipative force is large,
which cools down the system. When the system is cold, the random force term heats up the system. The thermostat
works well with equilibrium system.
But when come to the nonequilibrium system, i.e., Couette flow, because of the existence of the linear streaming velocity
gradient, when computing the dissipative force, the relative velocity may include a contribution from the shear field not just
the thermal velocity. When the shear rate is small, this might no be a big problem. But if I plan to apply an relative large shear
field, this contribution may take over the thermal part. I think this mechanism is suspicious under this situation. I hope someone
familiar with DPD could make a comment and give me some suggestion on this particular issue.
Best
Frank
Since the DPD thermostat uses the relative velocity
of 2 nearby particles, I would think it would still work
OK. Though as you say there is a 2nd order gradient
effect. You could also try using compute temp/profile,
which will measure the actual velocity profile, in conjunction
with a per-particle thermostat like fix langevin. I.e.
the compute temp/profile should subtract out the
streaming bias and just thermostat the remaining
thermal velocity. You could compare that to the
DPD pairwise thermostat.
Steve
Thanks Steve. Your quick reply is very helpful. Actually my result shows DPD work relatively well even at a substantial shear rate except near the boundary where the Lees-Edwards periodic boundary conditions applies on. The velocity profiles show a strange sudden change in slope as well as the temperature deviated dramatically from the value measure in the domain far away from the boundary (the temperature near boundary is much higher).
As we all know, when Lees-Edwards periodic boundary conditions assigned to maintain shear flow, any atom travel across that boundary will be subjected to a velocity change shear rate time boxsize in that dimension (i.e. \gammaprd[3]). My naive understanding is that the dissipative force between a pair straddling that boundary with be very large since the relative velocity will be in the order of \gammaprd[3]. It seems not physical. Is there a way to restrict the DPD interaction only between atoms in the simulation cell by switching off the interaction between any pair straddling the Lees-Edwards boundary.
Thank you for you kind reply sincerely
Frank
You are correct that there is an issue at the boundary. I hadn't
thought about that. Probably the right thing to do there is not
turn it off, but to have the ghost atoms be assigned a velocity
that accounts for the change in shear rate across the boundary,
similar to how their positions are adjusted by the periodicity.
Then there would be no difference in the relative velocities across
the shearing boundary versus in the middle of the box.
I'll have to look into the velocity communication logic to
see if that is easy to do.
Steve
Your comment is very helpful. One concern is if the ghost atoms are assigned a velocity taking out the change across the boundary, will that influence other types of interactions relying on the relative velocity who need actually count in the hrateprd part. The idea of turning off dissipative and random forces and only keep conservative force came from this paper A. Chatterjee, Molecular Simulation, 33, 2007. It’s quick and dirty, but I am not sure whether this treatment will bring unexpected effect on the system.
BTW, I am also looking in the source code. But I am a little confusing about where the ghost atoms’ coor and vel adjusted by the periodicity. Since when I looking at for example the pair_lj_cut.cpp, after looping over neighbors of one atom, the distances are simply
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delzdelz;
I see no condition structure to judge whether it’s the minimum image, so I assume you do it somewhere else, maybe embedded in the coord and vel communication of ghost atoms. Would you kindly pointed out the place in the source code where you adjust the coord. I look at Comm->border(), but didn’t see any adjustment.
Best
Frank
The code for coord adjustment across PBC is in the atom_vec_*.cpp
files. Nothing is done for velocities, but it could be for ghost atoms
if shearing is deformation is taking place. I can't think of many
places in the code where velocities of ghost atoms are used. I think
it would be correct to always adjust the ghost velocities across
a PBC due to the deformation gradient tensor.
Steve