Refactor bonus struct for ellipsoid atom style

[jibril-b-coulibaly]

Summary

I was looking into the ellipsoid atom style with the goal to implement granular contact between ellipsoids. I would like to better understand the rationale for the current bonus structure, which constains shape and quat, and discuss possible refactoring.

Detailed Description

The documentation mentions:

For the ellipsoid style, the particles are ellipsoids and each stores a flag which indicates whether it is a finite-size ellipsoid or a point particle. If it is an ellipsoid, it also stores a shape vector with the 3 diameters of the ellipsoid and a quaternion 4-vector with its orientation.

What are applications in which there is only a small number of ellipsoids and a large number of point particles, and for which defining the shape and quaternion only for actual ellipsoids saves up significant memory?

Atom style sphere does a similar thing where

If the diameter > 0.0, the particle is a finite-size sphere. If the diameter = 0.0, it is a point particle.

However, per-atom data that is not valid for point particles, e.g., omega is readily accessible and not located inside a bonus structure. What explains such difference?

Given that a quat pointer already exists in the Atom class, a shape pointer could similarly be defined for finite-sized particles and used for ellipsoids, thereby getting rid of the bonus struct and the complexity associated with handling the information (e.g., special checks and communication).

I would be happy to discuss more about the pros and cons of the current bonus structure vs having quaternion and shape for all ellipsoid atoms.

Thank you!

Further Information, Files, and Links

N/A

[stanmoore1]

What are applications in which there is only a small number of ellipsoids and a large number of point particles, and for which defining the shape and quaternion only for actual ellipsoids saves up significant memory?

A typical example would be ellipsoids in a stochastic rotation dynamics (SRD) solvent, see https://github.com/lammps/lammps/blob/develop/examples/ASPHERE/ellipsoid/in.ellipsoid and https://docs.lammps.org/fix_srd.html

[akohlmey]

The documentation is a bit misleading, since the ASPHERE package pair style for Gay-Berne describes an ellipsoidal generalization of the LJ potential. Also, there are multi-point "tri" and "line" particles that interact with LJ like groups of LJ point particles.

Please note that the bonus struct pointer is also used for atom style body and the pair styles in the BODY package.

[jibril-b-coulibaly]

Thank you very much for the prompt answers and clarifications.

I had found out that line, tri, and body also have bonus structures but their content seemed unique enough that it does not look like they could be consolidated into the default Atom pointers without significant changes. I would not wish to change the bonus structure for these styles.

I was curious about ellipsoid specifically because they only differ from spheres by the presence of the quat and shape pointer, and I thought these may be consolidated. I must clarify that I am looking at this from a "GRANULAR / DEM" perspective, in which case, there is usually no point particles: all particles are finite-size. Having shape and radius for any non-spherical finite-size particles could help define a circumscribed radius and a bounding box to help with contact detection.

I understand the memory concern for large systems with few actual ellipsoids like SRD. Given that bonus seems to be allocated and grown in chunks of 8192, any proc with more than 8192 atoms could potentially be wasting memory if the quat and shape pointers were defined for every atom and fewer than 8192 actual ellipsoids are defined. In the example for ellipsoid, there is 100 ellipsoids for ~15000 total atoms so that is a good illustration of the issue, thanks.

I wanted to clarify that aspect (and potentially propose a refactor) before implementing contact detection. I will do so with the current bonus structure.

Thank you

[jibril-b-coulibaly]

I would have one more question regarding the rationale for deciding which variable goes into a bonus struct and which should be a general Atom pointer.

For example:

• angmom
• omega
• radius
• torque

have a general definition and pointer in Atom, even though they are not relevant to point particles. These quantities are only used/defined for either tri, line, ellipsoid sphere or body. For those atom styles that have a bonus struct, this definition would cause unnecessary memory usage (similar to shape and quat for ellipsoids) if there are only a few tri, line, or ellipsoid among many point particles.

What are the main reasons why these particular variable are not part of a bonus struct?

[akohlmey]

What are the main reasons why these particular variable are not part of a bonus struct?

You have to understand that LAMMPS didn't happen instantly but has grown with contributions from different developers (some of which had quite different ideas of what is "obvious" programming). The attempts to remain consistent and even go so far to enforce consistent style and implementations are rather new (last few years and with extra force added during the pandemic, when there was more time to start more deeply reaching refactoring work).

That said, most of the properties you list are required by atom style sphere. The question whether to use bonus or add a new atom style is a bit of a personal choice probably augmented by the suggestions of @sjplimp, althrough I don't recall who actually came up with adding the bonus field and the different structs. It may have been someone else.
Please also note that we have now an alternate approach to add named per-atom properties with fix property/atom, that makes adding atom styles unnecessary in several cases.

When using atom style sphere, e.g. for DEM, it is not as likely to have a mix of point particles with extended particles. Also, atom style sphere predates the bonus structs.

[sjplimp]

The motivation for the bonus data structs for all the several atom styles (ellipsoid, line, tri, body) was to make
them efficient memory-wise when only a small fraction of the "atoms" in the system are ellipsoids, etc.
But a very large number of simple solvent particles are present. Which, for example, pair gayberne and other
colloidal pair styles support. See the blue Note on the pair gayberne doc page. The bonus data structs
then offer a large per-atom memory savings. And they are not particularly inefficient even for systems which
are all ellipsoids, etc.

Omega, torque, radius are all used by granular systems, which are typically all-granular, i.e. use of atom-style sphere.
So no motivation to make them part of a Bonus data struct. I can't recall which atom styles use angmom.

[akohlmey]

I can't recall which atom styles use angmom.

$ grep -l angmom src/atom_vec_*.cpp src/*/atom_vec_*.cpp
src/atom_vec_body.cpp
src/atom_vec_ellipsoid.cpp
src/atom_vec_tri.cpp

[jibril-b-coulibaly]

Thank you very much for these clarifications.

I will use the data structure as is, and interface it with pair_granular so that ellipsoids can be used for granular simulations using atom style ellipsoid.

If I may ask for advice, a circumscribed radius to the ellipsoid (i.e. the max of shape[3]) will be needed to facilitate contact detection: would you suggest creating a new bonus attribute, or using the radius pointer ? In the latter case, should theradius_flag be explicit set to false to indicate radius is only used for contact detection purposes in granular simulations?

Alternatively, granular contact detection between ellipsoids (x^2 + y^2 + z^2) and so-called super-ellipsoids (x^N + y^N + z^N) is nearly identical, I may create a new superellipsoid atom style geared towards granular simulations. This would encompass the ellipsoid capabilities but leave the current ellipsoid implementation unchanged.

I truly appreciate your guidance

[akohlmey]

If I may ask for advice, a circumscribed radius to the ellipsoid (i.e. the max of shape[3]) will be needed to facilitate contact detection: would you suggest creating a new bonus attribute, or using the radius pointer ? In the latter case, should theradius_flag be explicit set to false to indicate radius is only used for contact detection purposes in granular simulations?

I think using the radius pointer for a property that is not "the radius" is not a good idea. I would suggest to append a new property "maxrad" to the bonus struct for ellipsoids.

As for the super-ellipsoids, why not place a flag into the bonus struct for that and handle both in the same atom style? That we it may be possible to have both kinds of particles in the same simulation.

[jibril-b-coulibaly]

We could handle both of them in the same ellipsoid style.

That would also add 3 exponents to the struct (I wrote x^N + y^N + z^N but it's actually x^N1 + y^N2 + z^N3) that could default to 2, with an "all exponents are 2" flag.

This might be a bit invasive to check and change for existing code that assumes exponents of 2, e.g., to compute moments of inertia, or in the data file parser . But if that is preferred over a new atom style I'll take that route.

Thank you

[jtclemm]

I think using the radius pointer for a property that is not "the radius" is not a good idea. I would suggest to append a new property "maxrad" to the bonus struct for ellipsoids.

Although it would be a misnomer, one advantage of using atom->radius is that it already works with the current neighbor list classes (if the goal is to use granular pair styles).

[akohlmey]

Although it would be a misnomer, one advantage of using atom->radius is that it already works with the current neighbor list classes (if the goal is to use granular pair styles).

That is a valid point.

[jibril-b-coulibaly]

So far in #4008 , I have extended the ellipsoid data structure and properties (volume, inertia) so that the atom style ellipsoid can cover both regular ellipsoids and super-ellipsoids. I will now address contact detection and the handling of the radius.

Although it would be a misnomer, one advantage of using atom->radius is that it already works with the current neighbor list classes (if the goal is to use granular pair styles).

That is a valid point.

@akohlmey , based on @jtclemm 's comment, would you approve of using the atom->radius pointer to store the circumscribed radius for granular contact detection? With appropriate caution that it is not actually a radius (might echo some questions of #4054 with regards to flagging).

I currently implemented the circumscribed radius as a bonus quantity per your suggestion: what would be the alternative way to build the neighbor list from that information?

Thank you

[akohlmey]

@jibril-b-coulibaly I don't want to be the only person that has a say in this. At this point it would be good to have input from @sjplimp.

In fact, @sjplimp and I had a long discussion about refactoring the handling of per-atom data to make it simpler, future-proof and generally change the approach as we have a growing number of atom styles with overlapping feature sets.

I think the major question you need to ask yourself when re-using a per-atom attribute like "radius" in this case is: would it be conceivable (even if very improbable) that somebody might want to set up a simulation using a hybrid atom style that would include the same attribute and would there be a conflict? E.g. in this case, would there be a situation where one would want to use particles requiring either atom styles ellipsoid or sphere.

I don't think there is currently a way to build neighbor lists based on "bonus" information.

[jibril-b-coulibaly]

Thank you for outlining the major question, it helps me clarify the thought process and I will use it now and in the future.

I would argue in favor of using atom->radius for the following reasons:

1. the ellipsoid atom style generalizes sphere. I don't conceive a user using a hybrid of these two styles. This would give both omega and angmom attributes, and calculation and integration of those would likely be inconsistent and give inadequate values. This would also give radius and shape which similarly appear complicated to manage and not very useful together.

2. bpm/sphere is more or less a sphere with a quaternion. Likewise, I don't see a case for using a hybrid of this style with ellipsoid given the specialized use of bpm/sphere. bpm/sphere has an atom->quat attribute while ellipsoid stores its quaternion in the bonus strucut, which should also cause conflict / inconsistency.

3. Some atom styles that do not technically have a radius already use the radius flag and attribute: this is the case for line, tri and body, which use it as a circumscribed radius, similar to what we would want to do for ellipsoid.

   • Although this shows it is possible and is already implemented, it is unclear to me what it actually does.