Refactor bonus struct for ellipsoid atom style
jibril-b-coulibaly opened this issue · 15 comments
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
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
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.
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
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?
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.
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.
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
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
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 newbonus
attribute, or using theradius
pointer ? In the latter case, should theradius_flag
be explicit set to false to indicateradius
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.
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
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).
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.
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
@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.
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:
-
the
ellipsoid
atom style generalizessphere
. I don't conceive a user using a hybrid of these two styles. This would give bothomega
andangmom
attributes, and calculation and integration of those would likely be inconsistent and give inadequate values. This would also giveradius
andshape
which similarly appear complicated to manage and not very useful together. -
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 withellipsoid
given the specialized use ofbpm/sphere
.bpm/sphere
has anatom->quat
attribute whileellipsoid
stores its quaternion in thebonus
strucut, which should also cause conflict / inconsistency. -
Some atom styles that do not technically have a radius already use the
radius
flag and attribute: this is the case forline
,tri
andbody
, which use it as a circumscribed radius, similar to what we would want to do forellipsoid
.- Although this shows it is possible and is already implemented, it is unclear to me what it actually does.