Mueller Plathe Flow Updater

[Yuchen Xi]

Hi all

I am constructing a 3D system of rigid bodies, and trying to simulate a shear flow using the Mueller Plathe Flow. Below is the minimum script associated with the shear flow:

const_flow = hoomd.variant.linear_interp([(0,0),(5e3, 1*5e3)])

hoomd.md.update.mueller_plathe_flow(all, const_flow, hoomd.md.update.mueller_plathe_flow.Y, hoomd.md.update.mueller_plathe_flow.X, 100)

The velocity gradient is in y direction and the shear flow is in x direction.

While I run the script, the following warnings continuously appears in the console:

*Warning*:  After 100 MuellerPlatheFlow could not achieve the target: 4511 only 339.222 could be achieved.

I am not sure what that means. Has anyone had the same kind of issue?

Thanks

Lawrence

[Michael Howard]

Hi Lawrence,

Could you give us a little more information on your system? Do you just have rigid bodies (or is there solvent), and if only rigid bodies, how many / at what concentration? Here are a couple thoughts, some of which may or may not be relevant for you.

The warning likely means that you have set too high of a target flow, and the system is not able to swap enough momentum to reach it. This could just be because of your parameterization (in which case, choose a smaller target / longer time to reach it), or it could be because the method is not working very well (e.g., because you only have rigid bodies and there are not very many of them).

In general, I'm not sure how well HOOMD's implementation of Müller-Plathe's method will work with rigid bodies. I would imagine you only want to transfer momentum between the "rigid centers" (not all particles).

Just FYI, we have an alternate implementation of this method that uses the more "traditional" parameters, where you fix the number of swaps and the target momentum for each swap. We've found this has been easier for us to dial in the right shear rates, especially if you have some estimate of the viscosity already. See the azplugins package:

https://github.com/mphowardlab/azplugins/blob/master/azplugins/flow.py#L535-L592

Regards,

Mike

[Yuchen Xi]

Hi Mike

Thanks for getting back to me!

My current system has only 4096 rigid ‘rods’ where the rigid center is at the bottom. The system has a packing fraction of 0.84.

I do only want to focus on the rigid center, and i think that what you said makes sense. Maybe I should choose a longer time to reach the target flow (say 5e5 time steps)?

Lawrence

[Michael Howard]

Hi Lawrence,

No problem. That's a pretty high packing fraction, so I think you should be OK applying reverse perturbation to the rods directly.

Yes, I would try only using the rigid center group (this is where all the linear momentum of the rigid body is). You might then also try increasing the time to reach that target momentum. I believe the quantity you are specifying is the integrated (sum) momentum exchange, so you might try estimating the total momentum transfer as the amount you want to transfer per step times the number of steps, trying not to be too aggressive.

Also, how big is your box? You want to make sure your slabs are big enough to contain a sufficient number of rods for swapping, and you should choose n_slabs accordingly. You also need to make sure that the distance between the momentum exchange regions is large compared to the size of the rod. Last, make sure you keep the box a cube, or taller in the shear gradient direction:

https://doi.org/10.1103/PhysRevFluids.4.043905

https://doi.org/10.1063/5.0010697

Last, I would carefully check that the rigid constraint algorithms still work correctly with reverse perturbation. Maybe dump out a configuration after you run the updater for a while, and check that the rods & their constituents have preserved their original shape.

Regards,

Mike

[Yuchen Xi]

Hi Mike

My box dimension is 71sigma*71sigma*30sigma. The rigid rods are 5sigma tall in the z direction. I am trying to shear it with a velocity gradient in y direction and shear flow in x direction. Is 70 slabs a reasonable amount?

Lawrence

--
You received this message because you are subscribed to a topic in the Google Groups "hoomd-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/hoomd-users/rjHUD-rvFB4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to hoomd-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/af44f30f-3986-441c-8777-62f1bd3b651fn%40googlegroups.com.

[Michael Howard]

Hi Lawrence,

71 slabs (slab width 1.0) would probably be reasonable, but you should estimate the number of centers-of-mass in each slab based on your rod number density, the slab width, and the box cross-section. If you find there aren't many centers per slab, you could try decreasing the number of slabs a bit (say, 35 or 50 slabs).

Regards,

Mike

[Yuchen Xi]

Hi Mike

I've checked your azplugin reverse_perturbation guidelines, and there is a warning saying that it works only with NVE integration. Does that mean it is not compatible with hoomd brownian integration?

Lawrence

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/bd4efd74-4ef1-4299-9707-f58f287f9716n%40googlegroups.com.

[Michael Howard]

That's correct, and the implementation in hoomd is also not compatible with that integrator. Brownian dynamics is overdamped and has no momentum (the velocities are computed instantaneously from the forces), so the momentum-swap method cannot do anything. Other thermostatted integrators also will not work with shear flow because the temperature is computed assuming the net velocity is zero everywhere.

A main advantage of Müller-Plathe's method is that it is energy-conserving and hence you do not need to run it with a thermostat (NVE integrator works). You will potentially get temperature gradients, though, see Fig. 4 of the original paper:

https://doi.org/10.1103/PhysRevE.59.4894

Regards,

Mike

[Yuchen Xi]

Hi Mike

I kinda get it. But if I want to use Mueller-Plathe's method and dump a gsd file to visualize the motion of my system, which integrator in hoomd should I choose? 

Lawrence

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/85be2c92-08ba-43bc-81fc-7795c84f0b89n%40googlegroups.com.

[antoni...@googlemail.com]

Hi Lawrence,

Use the NVE integrator, for example: 

all = hoomd.group.all()
hoomd.dump.gsd(filename="trajectory.gsd",overwrite=True, period=1e2, group=all)
nve = hoomd.md.integrate.nve(group = all)
f = azplugins.flow.reverse_perturbation(group=all, width=1, Nswap=1, period=10, target_momentum=0.5)

Here, all particles are being integrated and swapped. Adjust this to your needs. Also, the reverse perturbation parameters might not be the correct ones for your system. Regardless if you use the azplugins implementation or the hoomd.md.update.mueller_plathe_flow you should be using the NVE integrator. 

Let us know if you have any more questions,

Antonia

[Yuchen Xi]

Hi Antonia

Thanks for your reply! 

After I use NVE integration, the Mueller-Plathe method and try to visualize the system, it is not "moving" at all. Below is the part of my script:

hoomd.md.integrate.mode_standard(dt=1e-6)
hoomd.md.integrate.nve(group=rigid)

const_flow = hoomd.variant.linear_interp([(0,0),(5e6, 0.1*5e6)])

hoomd.md.update.mueller_plathe_flow(rigid, const_flow, hoomd.md.update.mueller_plathe_flow.Y, hoomd.md.update.mueller_plathe_flow.X, 35)

hoomd.dump.gsd('gamma=0.1.gsd', period=5e6, group=hoomd.group.all(), overwrite=True)
hoomd.run(5e6)

I still continuously receive the following warning:

*Warning*:  After 100 MuellerPlatheFlow could not achieve the target: 29 only 0.00594962 could be achieved.

I am not sure whether the method is not suitable for my system or not.

Lawrence

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/839b70a2-91cf-403b-86bd-118f71d9b9f0n%40googlegroups.com.

[Michael Howard]

Hi Lawrence,

1. Does your system look like it's moving when you don't have the flow updater on? You have an extraordinarily small time step, and during the run, you only advance the system by 5 tau. You might want to see if you can increase dt.

2. You are dumping with a period equal to the total run time, so you are probably only writing the initial configuration. Trying writing more frequently or running longer.

3. Have you given the system an initial set of velocities (i.e., what is its temperature)? You wouldn't have needed to do this if you were using the BD integrator, but NVE requires that the system is initially equilibrated and thermalized.

I can't give advice on whether the parameters you are choosing or method are suitable for your system. I would again recommend installing and trying the version in azplugins, as I think those parameters are much easier to interpret (see Antonia's example) compared to the original paper describing the algorithm. That code will not generate the warning you're getting here, but you will still need to play with the parameters to see what generates the flow rate you want.

Regards,

Mike

[Yuchen Xi]

Thanks for your suggestion Mike! I will go and check mine.

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/8c27c762-5e8e-410f-8d54-12955adaa487n%40googlegroups.com.