SLLOD Algorithm

[astatt]

To implement the SLLOD algorithm for a homogeneous linear planar (Couette) shear flow, we need:

1. a triclinic box deforming from -45 deg to +45 deg and flipping back in a periodic way
2. each particle leaving the box in gradient direction at the bottom needs to get +shear_rate to its velocity, the ones leaving at the top need to get -shear_rate
3. A thermostat which applies a gradient field consistent with box deformation and thermal fluctuations relative to the gradient
4. And of course correct equations of motion.

Since these things are all tightly coupled together, only one class should be written: a SLLOD IntegrationMethodTwoStep. See TwoStepLangevinFlow for an example. The box can be updated to the current deformation based on the time step internally in the integrator, probably at the end of step 1, with a call to pdata->setGlobalBox(). When the integration of the equations of motion (i.e., update the positions) is done in step 1, the velocity wrapping for crossing pbcs can also be applied inside the integrator. If the box gets updated in step 1, then the forces will be computed before step 2 using the deformed box shape, correctly capturing the relative motion of the boxes.

The common convention in azplugins is to use z as gradient direction and x as flow direction (see RNEMD implementation).

The lammps implementation can be used as a model, see documentation and corresponding source code

[mphoward]

This is a good idea and would be very useful. A few comments:

1. a triclinic box deforming from -45 deg to +45 deg and flipping back in a periodic way

This will require a bit of thought because it requires remapping the domain decomposition. It probably means use of a snapshot in the C++ code, maybe as an updater so that it triggers first thing in the integration step. You would pull all of the particles back onto one rank, change the box, then broadcast them back out. I think there might be a method that does somewhere like this in the HPMC module, so we could try to look for that.

2. each particle leaving the box in gradient direction at the bottom needs to get +shear_rate to its velocity, the ones leaving at the top need to get -shear_rate

What happens to pair potentials like DPD that need PBC wrapping of the velocities? With this proposed implementation, I think there is no way to deal with this. Would we just document that this use case is not supported?

3. A thermostat which applies a gradient field consistent with box deformation and thermal fluctuations relative to the gradient

A profile-unbiased thermostat should be easy enough to add, as in the MPCD code. I'm not sure if there is a good way to inject this into HOOMD's core thermo computes, or if we just need to use an extra one for this integrator. DPD thermostat would also work, except for the reason listed above.

The common convention in azplugins is to use z as gradient direction and x as flow direction (see RNEMD implementation).

In HOOMD 3.0, I would like to rotate the coordinate system so that y is the gradient and x is the flow direction. This would allow methods to continue to function even in 2d systems (see this issue). Depending on when you work on this, though, I am OK if you want to use z for now since it shouldn't be too hard to tweak later.

[astatt]

You are correct, I have not considered pairwise thermostatting. See this lammps post: https://sourceforge.net/p/lammps/mailman/message/36415453/. After a quick search it looks like people generally ignore this issue.

We could think about using Lees-Edwards boundary conditions instead, since then the positions are modified for crossing pbc instead of velocities? If I remember correctly, we had said in previous discussions that LE boundary conditions would be very difficult to implement in hoomd.

[mphoward]

I believe that in Lees-Edwards boundary conditions, the velocities need to be modified when wrapped through the pbc. See these notes from David Kofke:

http://www.eng.buffalo.edu/~kofke/ce530/Lectures/Lecture24/sld024.htm

By doing the box deformation and wrapping the particle velocities when they cross the boundary, you should be doing this correctly.

My understanding is that the SLLOD part of this shear algorithm is actually only the equations of motion (and not the thermostat), so you could in principle run "SLLOD-NVE" not "SLLOD-NVT". SLLOD-NVE doesn't work in practice because you get viscous heating, but the DPD thermostat should be able to fix this just as well as a profile-unbiased thermostat. I think what Axel is suggesting in that post is that there is no way to disable the profile-unbiased thermostat in LAMMPS (it is hard-baked into the equations of motion), so using a DPD thermostat on top would be a bit weird... having 2 thermostats at the same time.

Unfortunately, pairwise thermostatting won't work with this implementation because you need to be able to compute the pairwise difference in the velocities of two particles. Like position, there is a minimum image convention for doing this in Lees-Edwards boundary conditions (ex: particles at the bottom of the box are moving left, but their images in the next box are moving right). However, hoomd's PairPotential just uses the BoxDim to compute dr and passes dv directly. Hence, I think there is likely no way to support this without modifying either the hoomd BoxDim directly, or by implementing a new version of the DPD potential that does support shear.

Given that this is specific to DPD, I'm OK with explicitly documenting that this is not supported in the plugin. If the code works well, we could always get it into HOOMD directly and fix it. There are very few parts of the code that might use the image of the velocity / velocity difference (DPD and MPCD being important ones for us though), so it wouldn't be too bad to change.

[astatt]

LAMMPS Pseudocode/algorithm

Lammps is using a Nose-Hoover thermostat, which is also updated/applied in two half steps. Tilt in xy direction, other tilts =0.

1. first half step of Verlet
   1. perform half step of Nose-Hoover chain thermostat variables update
   2. calculate SLLOD correction vdelu = H_rate*H_inv*v, this simplifies to vdelu_x = gamma_dot/Ly*v_y, vdelu_y=vdelu_z=0
   3. remove the "bias", e.g. the flow field, from the velocities
   4. apply the velocity scaling from the thermostat
   5. apply the SLLOD correction: v_x = v_x -1/2 \Delta t *vdelu
   6. restore the "bias" back in
   7. half step update velocities v = v + 1/2 \Delta t * force
   8. full step update positions x = x + \Delta t * v
   9. wrap the positions pbc
   10. if particle crossed pbc in y, apply velocity correction v = v +/- H_rate=gamma_dot
2. second half step of Verlet
   1. half step update velocities v = v + 1/2 \Delta t * force
   2. perform half step of Nose-Hoover chain thermostat variables update
   3. calculate SLLOD correction vdelu = H_rate*H_inv*v, this simplifies to vdelu_x = gamma_dot/Ly*v_y, vdelu_y=vdelu_z=0
      4. remove the "bias", e.g. the flow field, from the velocities
      5. apply the velocity scaling from the thermostat
      6. apply the SLLOD correction: v_x = v_x -1/2 \Delta t *vdelu
      7. restore the "bias" back in

[mphoward]

I talked with someone who is implementing some version of this. They ended up needing to modify core HOOMD to do it properly, so I vote that we close this as wontfix for now. We can reopen if someone wants it in future, and we think there's a way to do it in azplugins rather than HOOMD.

[astatt]

Interesting! Can you tell me who is working on it, I'd be happy to talk to them directly. I am not quite sure what the reason for needing hoomd core would be, so I want to talk to them. I have some code that is almost working I believe (minus DPD, MPI).

[mphoward]

It is one of Jeremy's former students (Deepak Mangal) who is now a postdoc with Safa Jamali. They are doing DPD with MPI. :-) I think they were only doing the Lees–Edwards part and not SLLOD.

[mphoward]

Jeremy & I will be implementing this in HOOMD proper, so I would suggest we not pursue it here & close this issue. We would welcome your help in adapting the code you have to HOOMD though!

[astatt]

Absolutely. What is the timeline on the hoomd implementation? Any ideas/plans?

One of my students is currently using this. As discussed above, right now the implementation does NOT work for pairwise thermostats (DPD) and also does not work for any kind of domain decompositions etc, so single GPU only. Not sure what of these shortcomings need to be addressed to go into hoomd proper.

[mphoward]

Awesome! We are currently porting the reverse perturbation shear flow code, so this would not be started until after that. We are also going to rework the MPCD cell list very soon to allow the box to deform, which is needed for this to work. The first step after that is to modify the BoxDim to be aware of the deformation rate. I think that will help fix most of the major issues you’re describing above!