freud
freud copied to clipboard
ClusterProperties should allow users to specify masses
Is your feature request related to a problem? Please describe.
Users who deal with molecules and wish to calculate radius of gyration or centers of mass with respect to periodic boundaries can use the ClusterProperties
class. However, this class's compute
method does not allow users to specify masses.
Describe the solution you'd like
We already have code for calculating weighted periodic centers of mass but this should be extended to the ClusterProperties
class so that radii of gyration and gyration tensor can also be computed correctly.
We should also add an example notebook that shows how ClusterProperties
can be used with molecules. The use of ClusterProperties
for this task isn't particularly obvious and we often have to explicitly show that feature to users who are asking for that calculation. We considered renaming ClusterProperties
after this confusion in #407, but ultimately decided against it. However, we can definitely improve the examples so that this capability is explicitly shown.
Additional context
PRs #553 and #677 can be used for references for the box's center of mass feature and demonstrate how to accept an optional array of length num_points
that has a default value of 1 for all particles.
As the gyration tensor Wikipedia article notes, the proposed mass-dependent quantity should be called the moment of inertia tensor and not the gyration tensor.
Although they have different units, the gyration tensor is related to the moment of inertia tensor. The key difference is that the particle positions are weighted by mass in the inertia tensor, whereas the gyration tensor depends only on the particle positions; mass plays no role in defining the gyration tensor.
Unfortunately, providing both computed properties for gyration/moment of inertia tensors leaves some ambiguity about how to handle the centers
property. We currently document that centers
represents the centers of mass -- but in fact there is no mass being used in the computation (mass = 1). Importantly, if this API accepts masses, we cannot use the mass-weighted cluster centers when computing the gyration tensor. The gyration tensor is computed relative to the "mean position," rather than the actual center of mass.
To simplify this a bit, I would suggest deprecating the gyration tensor property and change the entire ClusterProperties
class to use masses in all computations (defaulting to a mass of 1 for all particles). The gyration tensor is then simply the "default" behavior of the moment of inertia tensor when masses are not provided (and default to 1). The units are a little weird to document (would we include a factor of mass in the moment of inertia tensor?) but presumably users can figure out how to interpret the units if we document the computation itself.
Meanwhile, the radius of gyration Wikipedia article says:
Radius of gyration [is] defined as the radial distance to a point which would have a moment of inertia the same as the body's actual distribution of mass, if the total mass of the body were concentrated there.
In that article, the quantity clearly considers masses. That's another vote for changing from "gyration tensor" to "moment of inertia tensor" in my view. The radius of gyration would require an additional division by the sum of the clusters' masses before taking the square root:
R_g = np.sqrt(np.trace(self.moments_of_inertia, axis1=-2, axis2=-1) / np.sum(self.cluster_masses))
This aligns with the IUPAP definition on Wikipedia:
While we're at it, we should document some formulas (perhaps from Wikipedia articles) that match our implementations. ClusterProperties
is a fairly well-used feature of freud but it is lacking in mathematical documentation.
Is it necessary to change the Cluster
class or just ClusterProperties
, given that mass has nothing to do with actually determining the cluster network. Specifying masses to ClusterProperties
only will make it easier for users to try different mass vectors for the same cluster, without having to recalculate the network again. Though I realize that probably wouldn't be a common use case.
@AlainKadar Only ClusterProperties
would be affected. The classes are actually somewhat independent, if you already know the particle indices comprising each cluster. For example, in a polymer system where sets of particles are bonded, the indices of the particles in each polymer form each cluster.
Thanks @bdice. One other note is that the gyration tensor is normalized by number of points:
Whereas the inertia tensor is not normalized by masses:
The off diagonal elements are also negative in the inertia tensor (Comparison). I'm thinking its best to keep both attributes because of this.