Inducing flow in an SRD fluid

Dear LAMMPS users,

I am trying to simulate flow around dumbbell particles sandwiched between two plates and I have chosen an SRD fluid as a solvent. From what I understood from both the mailing list and the documentation, by using fix gravity one can introduce a pressure gradient and thus induce flow (this approach is suggested by a few articles, as well: 10.1039/C2SM26727F and 10.1103/PhysRevE.86.066703).

Since I am new to LAMMPS I decided to start simple by just adding the fix gravity to the in.srd.mixture example script provided by the software:

fix 1 big nve
fix 2 small srd 50 big 1.0 0.25 49894 lamda 0.1 &
radius 0.88 search 0.2 collision noslip shift yes 15135
fix 3 small gravity 0.5 vector 1 0 0
fix 4 all enforce2d

Yet, it seems as if nothing changes – if I run LAMMPS twice, once with this script and once with the 3rd fix commented out, I obtain the same numerical result. Am I missing something obvious?

Going through the mailing list for possible solutions, I came across this thread: which suggests one should modify the SRD fix to introduce the gravity term. Is this indeed true? I am running the 14th May 2016 version of LAMMPS under openSUSE 13.1 on a i3-3240 processor. Thank you very much for your time! Best regards, Rumen Georgiev MPI-Stuttgart


At least the streaming step for SRD has to be modified to include the ballistic
component of the forced flow. By reading the code in fix_srd.cpp it is not taken
into account.

Also, the most used boundary condition involves "ghost particles" in the wall,
this is also not part of the current fix.



Do any articles use gravity with SRD solvent? I don’t know

if anyone uses that kind of model. Why not apply gravity
to the large particles? Your script

is undoubtably complicated. I suggest starting simple,

e.g. pure SRD fluid or just big particles, see if you can equilibrate,

see if applying gravity does anything, see how bit the gravity

term needs to be, etc.


Do any articles use gravity with SRD solvent? I don't know
if anyone uses that kind of model.

It is fairly common, the original MPCD paper has a flow. Further studies by
Gompper (and collaborators) and the paper by Whitmer and Luitjen
are good references.

There are studies of particles in flows too (nanomotors by Tao and Kapral,
polymers by H. Stark and co-workers).




Yes, both articles I cited in the original mail ( and ) use this approach. One of them is even co-authored by you . I applied gravity to the big particles and indeed I observe motion in the direction of the acceleration. When I apply the fix to the SRD particles, however, nothing changes compared to the case where there is no gravity applied at all. Since I want a flowing solvent I think the pressure gradient (or its equivalent in this case, the gravity) should be applied to the solvent particles, not the big ones. I am able to equilibrate the system and I have tried applying the gravity fix to the SRD particles using several different orders of magnitude (from 0.01 to 10^4), yet to no avail – the velocity of the SRD particles does not change at all. I also tried using fix addforce on the SRD particles – the result was the same. I would be glad if you can give me some pointers. Mr. Pierre de Buyl has already suggested fix_srd.cpp needs editing for the gravity to be taken into account, but since you have a paper from 2012 about a similar problem I hope you can help me, as well. Best regards, Rumen Georgiev

The implementation of induced flow in that paper was

by the 1st author, Dan Bolintineanu (CCd). You can

contact him to see if he has a version of fix srd with

that option. That option was never released in main LAMMPS.

I doubt that he used fix gravity. As Pierre said, fix

srd handles the advection of the SRD particles itself,

so I think an added force term would need to be handled internal

to fix srd…



As Steve said, fix gravity won’t work for SRD, there are some minor changes to the code that we had to make for forced flow.

Attached is the version of fix srd that I used for those papers. Unfortunately there are MANY other changes/additions, some of which aren’t used or working, so it’s a bit messy. If you search for “FPOISS” in the code you’ll find the relevant parts. Basically, the SRD advection step is modified slightly to simulate the effect of a force acting on SRD particles (I believe only in the x-direction, but it’d be easy to modify for others).

I’m leaving for vacation in a few hours, so I’ll be out of touch, but if you have questions feel free to send them my way and I can try to help when I come back (6/21).

-Dan Bolintineanu

fix_srd_mod.tar (250 KB)

Dear Mr. Bolintineanu,

Thank you very much for the modified version of fix_srd. I have been trying different thing for the past month and I have managed to obtain a satisfactory result for a pure SRD fluid. There is still some slip, though – you can find a comparison between LAMMPS’ result and an analytical Poiseuille flow profile in the Velocity_profile.pdf attachment and a table of the system parameters I am using in the Parameter_table.pdf.

I essentially follow your approach described in the 2012 PRE paper and use a 40^3 simulation box, which is periodic in the x- and y-directions and has two walls in the z-direction. I thermalize my system for 300 000 timesteps and then take a snapshot every 10 for 100 000. Here’s also the fix_srd script that I am using:

fix 1 all srd 1 NULL 1.0 2.00 435435
mpc srd rescale no collision noslip reverse2 flow fpoiss 0.0005
thermo hecht 0.1 skipnone vp multi_poiss_all shift yes 7863467 inside warn

I have checked all the reverse rules and they result in more or less the same flow profile safe for some minor differences. Flow profiles obtained with reverse2 and reverse3 were actually identical, as far as I could tell. Therefore, the residual slip is not due to the bounce back rule, most probably. That led me to think the root of the problem might be the virtual particle rule.

The vp multi switch doesn’t work but you have already mentioned this is an issue as a comment in the fix_srd.cpp.

It looks as if multi_poiss_all and multi_poiss_partial yield the same result, which I guess corresponds to VP_multi (green circles) in Fig. 12 in your paper. So, here are my main questions:

  1. How can I properly couple the added virtual particles to the flow, that is, use a Gaussian with mean velocity corresponding to the Poiseuille profile at the particle’s position, rather than zero-centered Gaussian?
  2. Was this implemented in the fix_srd.cpp version you sent me?

I understand the physics behind it, but since my C++ skills are rather rudimentary I can’t tell if this feature has been added already or not and for the same reason I would rather avoid tempering with the code myself.

Thank you very much for your time!

Looking forward to hearing from you!

Best regards,

Rumen Georgiev

Velocity_profiles.pdf (111 KB)

Parameter_table.pdf (13.8 KB)


I think you’re right that one of the sources of the discrepancy is in the virtual particle rule, i.e. where you have to specify the virtual particle velocity to be consistent with a Poiseuille flow profile. This is implemented but commented out in the version I sent you – see lines 5210-5225 in fix_srd.cpp. I think it will work if you just comment that back in and recompile, but haven’t tried it. Note that it’s specific to the flow geometry (flow in x, walls at z = 0 and z = L).

That correction is generally impractical, since it requires you to know the flow profile a priori. In the PRE paper we did it just to show what it would take to fully eliminate the slip. Generally it’s a very small effect though. Looking at your plots of the velocity profile, I’m guessing there may be another source for the discrepancy (though I could be wrong). One useful check may be computing the actual viscosity of your SRD fluid in a simulation of a pure fluid in a 3D periodic domain, e.g. using the Mueller-Plathe algorithm in fix viscosity, at different shear rates. It should match the calculated value for a range of shear rates that includes the shear rate corresponding to your forced flow. An example of this is in Petersen, Matt K., et al. “Mesoscale hydrodynamics via stochastic rotation dynamics: Comparison with Lennard-Jones fluid.” The Journal of chemical physics 132.17 (2010): 174106.

Also, I’d recommend longer runs and possibly slightly larger velocity histogram bins, just to get smoother velocity profiles. Hopefully that will make it more clear whether there is slip at the walls or a difference in the actual and expected viscosities.


Dear Mr. Bolintineanu,

I did as you had suggested – uncommenting the lines in the source file and recompiling works: with velocity rescaling turned on I obtain a smooth parabolic curve which is fitted well to the analytical Poiseuille profile. Thank you!

When that was sorted out I moved on to the next part of my task – inserting a colloid particle in the flowing liquid. Since you have mentioned in another thread () that a work by Gompper (PRE 90 033314) seemingly resolves the issue of the long time diffusion coefficient you had been struggling with, I decided to try this approach but encountered a series of problems:

log.lammps (2.09 KB)

in.srd.colloid (986 Bytes)

data.srd.colloid (121 KB)


Looking at the Gompper paper, I’m not completely clear on some of the implementation details. But their results make it look like it solves most of the issues I was having with colloid diffusion in an SRD solvent, which is why I suggested in the earlier thread that it might be a good option in general. As you’ve noticed, it won’t work with the current fix srd in lammps, you’d have to do some additional code modifications. Certainly not a rewrite from scratch though.

Note that there are two algorithms in the Gompper paper, the angular momentum conserving one and the non-conserving one. If colloid rotation matters, then obviously you’d want the former. Even the translational diffusion coefficients I think will change if you don’t have angular momentum conservation. There are some important differences in how the fluid SRD particles interact with the colloid ‘mesh’ SRD particles between the two algorithms (see their paper).

I think the 3 modifications you list are certainly needed, but they imply some others, e.g. having the ability to distinguish between the different SRD particle types. Also, it’s not clear to me how the coupling between the colloid motion/rotation and the ‘mesh’ SRD particles is done, i.e. at what point are the positions of the colloid SRD particles updated based on the colloid center of mass motion, or how forces/torques on the colloid are modified in response to the SRD collisions.

Not to change the subject, but there’s a more recent paper by the same group (PRE 90, 032604) in which they abandon that approach altogether and stick to a more ‘traditional’ sphere collision approach, with what looks to be more detail on the implementation. As you can see from the paper, it’s still a fairly complicated story, but this does seem more amenable to the current LAMMPS implementation.

Let me know what you think of all these different approaches. I haven’t looked at the details of these and related recent papers, but it seems like it’s a fairly active topic in the last couple of years (at least in the very small world of SRD/colloid coupling). I’d be interested in modifications to fix srd in LAMMPS that have some of these options, and may be able to help a little, let me know what you decide.


Dan Bolintineanu