Help on basic usage example to calculate shear viscosity according to the mueller_plathe_flow updater

[Rasmus "Termo" Lundsgaard]

Dear all.

I am a new hoomd user (but long time MD user), and I have come to hoomd to test out the possibility to "measure" the shear viscosity with the mueller_plathe_flow updater.

I am a bit confused on how to use the updater and how to measure the viscosity results from the algorithm.

After setting up the updater, I can start the run (the updater will run for every timestep) and get the "summed_exchanged_momentum" after the run.

As I understand from the method, the viscosity is proportional to this summed_momentum_flux/time and the gradient chosen for the const_flow variant, right?

So how many timesteps is needed for the run with the updater? wouldn't it make sense to log the summed_momentum_flux by some desired period?

Can someone provide a minimal example of how to calculate the shear viscosity using the method?

Regard

Rasmus

[Michael Howard]

Hi Rasmus,

RNEMD simulations can be a bit finicky to setup, as you need to choose appropriate parameters for your system to ensure that you are in a regime where the desired steady flow can be generated and the system has linear response. In general, you would impose a momentum flux, measure the velocity profile in the system, and extract a shear rate from the velocity field. A plot of the momentum flux vs. shear rate for a few different fluxes should be linear, and its slope as the shear rate goes to zero gives the zero-shear viscosity. You need to run the simulation long enough for the steady flow field to develop and then for you to accurately measure everything you need; this will depend on your model. I would recommend reading a bit from Muller-Plathe's original publication that discusses some of these issues. We have also recently done a critical review of the method, which is up as a preprint--you might find Sec. II of that helpful.

I believe that the HOOMD RNEMD updater works a little bit differently than the standard RNEMD algorithm, as you specify a target flow and then it tries to adjust the RNEMD parameters for you. The author of that method (Ludwig Schneider) would need to comment for you on how that actually works and the design of it, as I prefer to use my own implementation that is closer to the original method.

Regards,

Mike

[Ludwig Schneider]

Hi Rasmus,

Mike is right. The time of your simulation needs to be long enough that
you get a stable velocity profile. The slope of this profile gives you
the strain of your system. Measuring the velocity profile takes some
extra effort. Maybe you can do it with a post-processing step.

I can add the comment about the target flow: the target flow is the
integrated input stress. So for a constant shear flow you choose a
linear increase and slope defines the stress in your system. Makes
afterwards sure that you this stress is always imposed in your system.
(If that is impossible, you get a warning. That can be the case if your
stress is too high.) The observable 'summed_exchanged_momentum' is the
actually exchanged momentum. If everything is okay it is identical with
the target flow. (There is a constant factor between those, which is the
area of the simulation box perpendicular to the shear gradient.)

The advantage of the target flow instead of the original algorithm is
that (i) the stress input fluctuates as little as possible and (ii) you
can do other shear experiments than constant stress e.g. oscillatory shear.

The ratio between stress and strain yields the viscosity.

Regards,

Ludwig

On 29.01.19 16:31, Michael Howard wrote:
> Hi Rasmus,
>
> RNEMD simulations can be a bit finicky to setup, as you need to choose
> appropriate parameters for your system to ensure that you are in a regime
> where the desired steady flow can be generated and the system has linear
> response. In general, you would impose a momentum flux, measure the
> velocity profile in the system, and extract a shear rate from the velocity
> field. A plot of the momentum flux vs. shear rate for a few different
> fluxes should be linear, and its slope as the shear rate goes to zero gives
> the zero-shear viscosity. You need to run the simulation long enough for
> the steady flow field to develop and then for you to accurately measure
> everything you need; this will depend on your model. I would recommend
> reading a bit from Muller-Plathe's original publication

> <https://doi.org/10.1103/PhysRevE.59.4894> that discusses some of these

> issues. We have also recently done a critical review of the method, which

> is up as a preprint <https://arxiv.org/abs/1811.04097>--you might find Sec.

[Rasmus "Termo" Lundsgaard]

Dear all

I now found some time to continue my test of using this algorithm to model the shear viscosity of different systems.

I am working on GC systems with the Mie potential. For single particle systems I can get the flow profile to develop fine with a system of single GC particles, but with bonded particles I can not get the flow profile to develop at all (see attached figures).

Is there something I am missing in the setup for systems with bonded particles? When swapping of velocities happen, should it maybe happen all particles in the molecule that has the lowest and highest value?

Best regards

Rasmus

[Ludwig Schneider]

Dear Rasmus,

great that you could get the MP flow working in your system.

By the look of your profile, I would say the system did not reach a
steady state, yet. The peaks in the active zones are typically less
pronounced after the system reached a steady state.

For your bonded system, I just want be sure that you do not form a
percolating network with you bonds, do you? A network is a solid, it is
not expected to flow. Also I hope your bonds don't make your system
glassy. Glasses flow on time scales we probably cannot simulate. Just
that we are on the same page.

If you don't have a network, but something like liquid, linear polymers,
stars, dendrimers or so, it is expected for the viscosity of such system
to be significantly higher compared to the simple liquid.

I would say, you do not see profile, yet, because the slope of your
profile is at least an order of magnitude smaller than in your liquid.
This small slope is probably hidden in the thermal fluctuations of your
system.

This is absolutely normal: the thermal velocity is larger than the
collective flow motion, I see it in my polymer simulations as well. You
can get around it, by (a) choosing a larger system to get some
self-averaging or (b) by time averaging in the steady state.

I would suggest that you let your simulations run longer and try to
average the profiles afterwards.

Hope, I could help.

Ludwig

> [image: 1_center.png] <about:invalid#zClosurez>

[Rasmus "Termo" Lundsgaard]

Not it si for now very simple liquids as coarse grained models of alkanes. They equilibrium densities are spot on, so I trust the simulation is correct as such. I can now see that my problem is related to liquid phase simulations. If I do MP for 1 segment compounds I get the nice flow profile for those that are in gas phase but not when going into liquid state. (see videos in this link https://www.dropbox.com/s/8gncpbdk4wq3cxg/MP-videos.zip?dl=0). this of course suggest as you say that simulations have to run longer. I try for now to redo the results for hexane from The_shear_viscosity_of_molecular_fluids_A_calculation_by_reverse_nonequilibrium_molecular_dynamics, where they state that running with timesteps of 2fs for 1ns they reached full flow profile. I have now run for 10ns and I still don't see a flow profile?? If I push the slope higher because maybe my flow-profile simply is to little to see in the thermal fluctuation even when time-averaging I get alot of:

*Warning*:  After 100 MuellerPlatheFlow could not achieve the target: 0.056088 only -1.50443 could be achieved.
*Warning*:  After 100 MuellerPlatheFlow could not achieve the target: 0.056094 only 3.04012 could be achieved.

Any ideas to why I can not see a flowprofile? As my gas phase simulation shows it seems that my code is correct for making the flowprofile and getting the data out at least.
 

[Ludwig Schneider]

Dear Rasmus,

I am not familiar with your specific system. So I cannot really help you
to figure out some of your particles are not flowing as you expect them
to. I work mostly with softer models. I don't know if you can expect the
correct dynamic properties if you match the densities. Also, which
thermostat (if any) are you using in combination with the MP_flow? Many
thermostats are incompatible with shear flow because they do not
conserve momentum locally. I use DPD as implemented in HOOMD, but I also
know that Lowe-Anderson (not implemented in HOOMD)  would also work.

However, I can explain you the warning:

The warning is triggered because the MP_flow cannot fulfill your stress
target. It performed 100 swap updates at each time step and still
couldn't reach the desired target.
The simulation is continued anyways but you have to expect that the
actual stress in your simulation is lower than you specified as target.
You can inspect the real stress via the "summed_exchanged_momentum"
quantity provided by the MP_flow class. If your stress is too high for
the system, the actual stress as logged would be lower than your target.
(Note that the cross-section area is a factor between your flow target
and "summed_exchanged_momentum".

The reason for this deviation occurs most of the time if the maximum
momentum in the max. slab is smaller than the minimum momentum in the
min. slab. Swapping these momenta than actually reduces the stress
instead of increasing it as desired. This is an inherent limitation of
the MuellerPlathe RNEMDS scheme. The easiest way to work around this, is
to increase the number of particles in the max. and min. slabs by
increasing the system size (or using less slabs for your system if
possible).

There is also the slight possibility that the MP class uses a too tight
tolerance between target and actual swapped momenta. The result would be
a back and forth swapping of momenta. The "summed_exchanged_momentum"
would fluctuate around your target flow in this scenario. In this case,
you can use the "getFlowEpsilon" and "setFlowEpsilon" in order to adjust
the tolerance to higher values. However, I never encountered this
problem before. It may also be related to very small system sizes.

I hope you are on the way on tracking down the problem in your simulations.

All the best,

Ludwig

On 29.03.19 09:22, Rasmus "Termo" Lundsgaard wrote:
> Not it si for now very simple liquids as coarse grained models of alkanes.
> They equilibrium densities are spot on, so I trust the simulation is
> correct as such. I can now see that my problem is related to liquid phase
> simulations. If I do MP for 1 segment compounds I get the nice flow profile
> for those that are in gas phase but not when going into liquid state. (see
> videos in this link
> https://www.dropbox.com/s/8gncpbdk4wq3cxg/MP-videos.zip?dl=0). this of
> course suggest as you say that simulations have to run longer. I try for
> now to redo the results for hexane from
> The_shear_viscosity_of_molecular_fluids_A_calculation_by_reverse_nonequilibrium_molecular_dynamics

> <https://www.researchgate.net/publication/27269895_The_shear_viscosity_of_molecular_fluids_A_calculation_by_reverse_nonequilibrium_molecular_dynamics>,

[Rasmus "Termo" Lundsgaard]

Dear Ludwig,

thank you for you time.

Actually it seems that in my case I consistently get a to high stress compared to target. here a small dump from my log, where you see the iter warning and then my dump of actual stress and the target stress at timestep.

*Warning*:  After 100 MuellerPlatheFlow could not achieve the target: 0.268826 only -1.58595 could be achieved.
*Warning*:  After 100 MuellerPlatheFlow could not achieve the target: 0.268827 only 3.56478 could be achieved.
stress   0.2788 0.2600
stress   0.2788 0.2700
*Warning*:  After 100 MuellerPlatheFlow could not achieve the target: 0.288803 only 3.88521 could be achieved.
*Warning*:  After 100 MuellerPlatheFlow could not achieve the target: 0.288804 only -1.58192 could be achieved.
stress   0.2880 0.2800

As it can be seen here if the new stress-target already is within flow_epsilon nothing is of course updated (the 0.2788 is within flow_epsilon for the new target of 0.27, as it did this overshoot in the prior timestep.), so if I increases flow_epsilon I just get alot of timesteps with no updates. reducing flow_epsilon leeds to alot of 100 iter warnings 2 questions arise from this:

• How do I work against this overshooting? I tried to reduce timesteps, but it seems to follow the same pattern of overshooting the stress induced.
• How can the iteration iterate past the target without stopping, maybe and iteration stop should be introduced if fabs(actual stress) > target ?

Best regards

Rasmus

[Rasmus "Termo" Lundsgaard]

Got it now

so Issue was that slope was to low and even 1st iteration overshoots (and then it continues until 100), which confused me with the 100 iteration warning in output. By loosing epsilon and setting the slope higher and disregarding the first inital timesteps with many warning because of the 1 iteration overshoot, it can reach a full developed flow profile in less than 1ns...

Thank yo for getting me in the right direction - I guess part of my problem is because of the coarse grained system and for that reason some rather large and heavy particles.

[Ludwig Schneider]

Dear Rasmus,

great that you could solve the issue for you.

I guess what happens is, that the algorithm tries to swap momenta back
and forth between the slaps in the first time step. But since it
overshoots the target in any case, doing 100 of these swaps does not
solve the issue, hence the warning. As long as your flow target is
reached at later times, you can safely ignore the warning.

As an alternative, if you can increase your cross section area A
perpendicular to the shear gradient, the impact of a single momentum
swap becomes smaller, relatively. This needs a bigger system and thus
more computation time of course.