glotzerlab / hoomd-blue

Molecular dynamics and Monte Carlo soft matter simulation on GPUs.
http://glotzerlab.engin.umich.edu/hoomd-blue
BSD 3-Clause "New" or "Revised" License
334 stars 130 forks source link

A more flexible box resize updater #1640

Closed rpcollanton closed 6 months ago

rpcollanton commented 11 months ago

Description

From a close read of the source for the BoxResizeUpdater, it seems to me that it resets the dimensions of the entire box every time it updates the box, even for dimensions that are not different between the two boxes. This is fine in an "NVT" ensemble calculation, but does not seem compatible with a mixed ensemble calculation such as allowing two dimensions of the box to be controlled by a barostat. With the current framework, is there a way to adjust a single dimension manually while allowing other dimensions to be managed by a barostat?

Additionally, as far as I can tell, there is no current support for adjusting box dimensions independently with a more complicated function than a linear interpolation between two endpoints. For example, uniaxial elongation at constant volume isn't supported, except for by stopping and starting the simulation. This could also likely be accomplished with a custom updater. However, in both cases, given the need to rescale every particle's coordinates, leveraging GPU resources would be ideal, and I don't think you can do that with a custom updater (I could be wrong).

Proposed solution

I am almost finished writing a solution for this through modification to the C++ layer of HOOMD-blue, and would be excited to discuss whether including this in HOOMD-blue would be appropriate, and how to tweak the design to be in accordance with HOOMD-blue's design philosophy.

The current form of my solution works as follows: the Python layer takes as input up to 6 variants, and dimensions to which these variants correspond. Any dimension not specified is untouched by the box updater (rather than maintained at its t=0 value). Hence, if something else (say a barostat) touches the box dimension, this isn't undone by the box updater. When triggered, the box updater sets each paramater value i to box0[i]*variant_i(timestep)/variant_i(0). The scale factor is normalized to the variant's value at time zero such that the box starts at its original size.

This supports varying any set of box dimensions with complete independence, but changes the philosophy from "interpolating between an initial and final box" to "scaling the original box".

To not break backwards compatibility, this could be implemented in a few different ways -- 1) as its own separate Python class with its own separate C++ class, 2) as a switch within the current "hoomd.update.BoxResize" Python class that defaults to the original behavior and accepts the original arguments, but drops into different C++ classes, 3) as two separate Python classes that use the same C++ class (the new one) in different ways (since the original behavior is a special case of the general behavior I am suggesting), or 4) as a single Python class and single C++ class, with the same switch as mentioned in (2).

I am not sure of a clean way to accomplish the "switch" between behavior that I am suggesting, but it seems that it would absolutely need to be there if a single Python class is to be used, to support backwards compatibility.

Ultimately, if we decide to proceed, I have a lot of the code written for this and it could be reorganized accordingly. Looking forward to any thoughts or feedback!

Additional context

No response

joaander commented 11 months ago

Note: 1) Custom updates can call hoomd.update.BoxResize.update to get the performance of BoxResize with the flexibility of custom code. 2) Variants need not be linear.

That being said, BoxResize is limited. I support adding a new class (at the Python and C++ levels) that enables new use-cases such as those you propose. Many optional "mode" options can be confusing to users. I would suggest refactoring to combine BoxResizeUpdater and BoxResizeUpdaterGPU into one class that dispatches the scaleAndWrapParticles method appropriately based on the ExecutionConfiguration. Then a class that implements new functionality for choosing the box can subclass BoxResize.

Based on your proposed method, I suggest naming it BoxScale. box0 should be a parameter that the user provides and not inferred from the first step the updater is triggered. This enables the method to function appropriately when jobs are continued over multiple invocations of python.

How would you express the example use-case of uniaxial elongation at constant volume with the proposed API?

joaander commented 11 months ago

On second thought, I don't see how scaling could possibly work for the box tilt factors. Say the user starts with a cubic box (all tilts 0). variant * initial_tilt will never produce any tilt except 0.

rpcollanton commented 11 months ago

I hadn't thought about that, you're absolutely right. I think there might be a good way to tweak my proposed solution. The main advantage is not necessarily interpreting the box updates as "scaling", but the ability to independently adjust each parameter with its own variant. I am going to work through a couple of ideas and get back to you.

joaander commented 11 months ago

HOOMD 2.x used this approach: https://hoomd-blue.readthedocs.io/en/v2.9.7/module-hoomd-update.html#hoomd.update.box_resize . Many users found it very confusing to create many different variants which is why we simplified it to 1 variant and interpolation starting in 3.0.

For cases like your example "uniaxial elongation at constant volume", I don't think the 6 separate variants solves the issue in a clean way. Those 6 variants would need to be updated carefully in a co-dependent way to achieve the target.

Ultimately, what we need is a more flexible BoxVariant system that describes the box parameters as a function of time. Subclasses of this could implement linear interpolation, uniaxial compression, and other use-cases as needed. Certain subclasses could set ignore flags to retain existing dimensions as you propose. In this design, I would not implement a new BoxScale class. Rather I would update the BoxResizeUpdater C++ class to the new API and implement a backwards compatible API in the Python class - deprecating the old API to be removed in a 5.0 release.

rpcollanton commented 11 months ago

I was thinking about the general idea of a "vector-valued variant". Would there be any sense in generalizing it like that, rather than as a specialized BoxVariant? Either way, I think this is a good and valuable path forwards.

How would the case of varying only a single parameter and leaving the other parameters untouched look? This is not the same as maintaining the other parameters as constant, since there may be a barostat adjusting them. Maybe this could be accomplished with a list of 6 flags given to the Python class that are passed into the C++ class, where each flag indicates whether to vary that box parameter, and the BoxResizeUpdater ignores the value of the BoxVariant for all the non-varied parameters?

I am interested in implementing this functionality once we decide on something.

joaander commented 11 months ago

Yes, we would welcome your contribution. We discussed the options for this internally. We would like to implement this with a general VectorVariant that could possibly be used elsewhere in the code (e.g. NPT stress tensor). VectorVariant would work like Variant, except it would return const std::vector<Scalar>& instead of Scalar. In this way, all the possible box moves can be naturally expressed in a variety of pre-built (or custom VectorVariant subclasses). The box vector would be the standard notation [Lx, Ly, Lz, xy, xz, yz]. Box may need a to_vector helper method to produce the arrays in that format.

To support setting only some of these quantities, BoxResize could take an array of 6 flags to control which it sets.

Thanks for being willing to work on this. Feel free to reach out with any questions you might have. You may find it useful to open a draft pull request so we can see your work in progress code and offer suggestions based on that.

Once a linearly interpolating VectorVariant is available, then we would rewrite the BoxResize C++ code to use it. For backwards compatibility, users can still provide box1, box2, and variant. At the Python level, this would internally create the linearly interpolating variant and trigger a deprecation warning. Users can instead provide the box_variant directly.

joaander commented 11 months ago

@rpcollanton do you plan to make these improvements?

rpcollanton commented 10 months ago

Yes, I would be happy to and I plan to!