Converging-diverging channel flow simulation

[Spela Ivekovic]

Hi,

I would like to run a fluid dynamics simulation with a converging-diverging channel, where two dimensions have periodic boundary conditions and one doesn't.

Are there any examples of code using custom geometries that I could look at for guidance? I'm looking for ideas on how to construct the channel, i.e., define the heights to make it converging-diverging, define channel length, throat position (where it changes from converging to diverging) and boundary properties.

Is this something that would be better done with a separate software tool and then imported into an OpenMM simulation in the initialisation phase?

Thank you,

Spela

[Peter Eastman]

Hi Spela,

Strictly speaking, OpenMM doesn't support systems that are periodic in some dimensions but not others. In practice, though, you can make it appear non-periodic in one dimension by simply making the periodic box sufficiently large in that dimension so that periodic copies never interact with each other.

The next question is how you want to describe the channel. Is it defined by a set of particles in fixed positions? Or some sort of analytic potential function? (The latter could be implemented with a CustomExternalForce.) Or something else?

Peter

[Spela Ivekovic]

Hi Peter,

the channel wall is defined by a set of particles in fixed positions. How can I specify a particle with a fixed position (one that won't move during the simulation)?

The other questions I have are (I apologise, this would have been better placed in the help section of the forum):

- I need to perform measurements on my channel. For example, I'd like to be able to do cumulative averaging as well as measure mass flux through the channel at a specific location. Since passing the positions and velocities of molecules to the CPU at every time step is time consuming, I presume some of this functionality has probably been implemented on the GPU and measurement functions already exist?

- Can you apply an external force to the equations of motion of the molecules (for example, a gravitational force) in order to create the flow through the channel?

- I'd like to recompile the OpenMM library from the source code so that I can customise some of the code for my simulation without having to request new features from the OpenMM team. Are there instructions available for how to compile the libraries from source?

Thank you very much.

Spela

[Peter Eastman]

There's no way to keep a particle absolutely fixed, but you can do a pretty good approximation of it. First, give it a very large mass (like 1000 or 10000) so forces on it will have very little effect. Then use a CustomExternalForce to apply a harmonic restraining force binding it to the desired position. The API docs for CustomExternalForce give example code for doing that.

> I need to perform measurements on my channel.

I'd start out by doing this in the obvious, straightforward way. Just download the coordinates every time step (or less often if that's ok) and compute whatever you need. There's some overhead from doing that, but as long as your system isn't too small, it won't be that bad. You can always do something more sophisticated later if it turns out to be necessary, but there's a good chance it won't be.

> Can you apply an external force to the equations of motion of the
> molecules (for example, a gravitational force) in order to create the
> flow through the channel?

Yes, using CustomExternalForce. A gravitational force in the z direction would be implemented by something like

CustomExternalForce("g*z")

where g gives the strength of the force.

> Are there instructions available for how to compile the libraries
> from source?

Yes, the manual has a whole chapter on it.

Peter

[Spela Ivekovic]

Hi Peter,

thank you very much for your answer.

I have successfully resolved all of the above.

For the benefit of anyone else reading this forum, I'd just like to report that I attempted using the CustomExternalForce to apply harmonic restraining force to the channel wall as well as to apply the gravitational force to the fluid molecules moving through the channel and have encountered two setbacks:
(1) CUDA platform does not allow more than one instance of the CustomExternalForce. This means you cannot define the harmonic restraining force as one instance and then gravitational force as another instance, they have to be combined together into one CustomExternalForce instance only.
(2) CUDA platform does not support more than 4 parameters in the energy equation. This means that I could not combine the two forces into one CustomExternalForce expression either because the harmonic force uses 4 parameters (k,x_0,y_0,z_0) and the gravitational force another one, say f_g.

In the end, this was resolved by not using the harmonic restraining force but instead increasing the masses of the wall molecules to make them stay put and only applying the gravitational force as a CustomExternalForce to the fluid molecules.

Spela

[Peter Eastman]

A simpler solution is to use OpenCL instead of CUDA. It supports multiple CustomExternalForces, and also supports an arbitrary number of parameters for each one.

Peter

[Charles Brooks]

This follow-up is regarding the limitations in the Custom Forces in OpenMM on CUDA. This particular post was regarding CustomExternalForce limitations. My questions are as follows: 1) Does this restriction still apply for CUDA, 2) does it extend to all Custom Forces, e.g., CustomNonbondedForce, CustomTorsionForce, CustomCompoundBondForce, etc., and 3) is it documented somewhere?

We have been extensively using the CustomForce apparatus to bring CHARMM functionality into the CHARMM/OpenMM interface and have recently noted problems (on CUDA) of an implementation of forces using the CustomCompoundBondForce, where the code seems to work on the Reference and OpenCL platform, but fails in somewhat different places depending on the compiler platform used to build the code linking to the OpenMM/5.2 pre-compiled libraries (GCC or Intel based compilers).

Thanks for any comments/clarifications.

Charles Brooks

[Peter Eastman]

This restriction no longer applies. Since version 5.0, the CUDA platform supports as many forces as you want of each class.

Peter

[Charles Brooks]

Thanks Peter. We have, then a curious issue we'll continue tracing. We have a case where we have implemented multiple CustomCompoundBond Forces, and as noted, they always seem to work on the OpenCL platform but we see varied results from the CUDA platform, depending on the compiler, as I noted, and also, apparently, which hardware installation. I'll let you know.