Universe.dimensions: enhancements, guarantees, automatic conversions

[zemanj]

Scope

This issue is a follow-up of #2190 discussing possible enhancements related to Universe.dimensions and, likewise, coordinates.base.Timestep.dimensions.

Questions to be addressed

• [ ] What is considered a valid box in MDAnalysis?
• [ ] How should missing box information be represented and how can it be distinguished from open systems?
• [ ] How should semi-periodic boxes be represented?
• [ ] (How) can we distinguish between semi-periodic open systems and semi-periodic closed systems (i.e., with walls)?
• [ ] Do we need support for true 1D or 2D systems?
• [ ] How can we detect the different cases?
• [ ] Should MDanalysis try to convert invalid boxes into valid ones automatically? If yes, how?
• [ ] How can we guarantee that round-tripping a trajectory through MDAnalysis doesn't change or degrade the box information in that trajectory?
• [ ] Is there a benefit of storing box information in double precision?
• [ ] If we have answers to the above questions, can we make guarantees about what Universe.dimensions returns?
  • [ ] How can we enforce such guarantees?
  • [ ] Is there a central point of enforcement for both Universe and Timestep?
  • [ ] Can we use such guarantees to simplify box checks in individual methods?
• [ ] Is there a central point in the docs where these things can be summarized?

Other related questions

• [ ] What is the impact on volume computation? Does it have to be adjusted?

If you have additional questions in mind, feel free to add more points to the list!

[richardjgowers]

So currently, we support triclinic boxes, and an orthogonal box is a special case of that which allows some optimisations. There's only one place dimensions is stored (like positions), all other access points are just aliases, ie [AtomGroup/Universe].dimensions == Timestep.dimensions

I think a nice theoretical solution to most of your questions is to go object-oriented, and make the box have methods which figure out periodic images, so an implementation of calc_bonds would look like:

def calc_bonds(pos1, pos2):
    vectors = pos1 - pos2
    minimum_vectors = box.apply_minimum_image_convention(vectors)

    return norm(minimum_vectors)

So here we're passing off responsibility to a Box object, they decide how that's done. Then if someone wants to do some weird semiperiodic/wall/whatever they have all the power to do that.

In terms of actual implementation... obviously its trickier as we work down at the pure C level. One option would be to rewrite c_distances in cython (annoying but possible...) then Box is a cython extension type, and all is well. Then our signature for the above function is:

def calc_bonds(float[:, :] positions1, float[:, :] positions2, Box box):

If Cython is slower (hopefully not) we could either keep a menagerie of functions (as currently) and have some sort of case select which looks at the type of the Box and selects the appropriate pure-C function (ie essentially not binding the function to the object) but that's ugly as sin. Could probably do some C++ binding instead, have C++ classes, but like came up elsewhere, this is a python project (ideally)

The idea of having .dimensions become an object could also lead to some nice API possibilities where people could be rolling their own distance code at near numpy-like performance as we've exposed the minimum image calculation neatly...

u = mda.Universe()

vecs = u.atoms[:10] - u.atoms[10:20]

min_vecs = u.dimensions.minimum_vectors(vecs)

tl;dr Cython?