Hexatic normalization in docs doesn't match implementation

[bdice]

@Plastikschuessel noted that the documentation for the Hexatic order parameter doesn't match the implementation. The docs say that the formula normalizes by 1/n, but the implemented code normalizes by 1/k. In many systems (particularly hexagonal and square crystals), k=n, but in some cases (e.g. 2D quasicrystals), the order is 12-fold even though any given particle has fewer than 12 neighbors.

freud/freud/order.pyx

Lines 250 to 251 in 3f44951

	:math:`\psi_k \left( i \right) = \frac{1}{n}
	\sum_j^n e^{k i \phi_{ij}}`

freud/cpp/order/HexaticTranslational.cc

Lines 38 to 45 in 3f44951

	if (m_weighted)
	{
	m_psi_array[i] /= std::complex<float>(total_weight);
	}
	else
	{
	m_psi_array[i] /= std::complex<float>(m_k);
	}

A temporary workaround to achieve the desired behavior is to enabled weighted=True in the constructor. In ball queries and nearest neighbor queries (but not Voronoi queries), the weight for each bond defaults to 1. If weighted=True, the normalization will divide by the number of neighbor bonds (a sum with 1 for each bond).

To resolve this issue, we should verify the literature to see which convention is more common and consider changing to normalization by 1/n instead of 1/k.

[bdice]

We should also clarify two things:

1. The Hexatic order parameter returns the complex valued OP. This is not rotationally invariant -- its phase has physical meaning (the orientation of the crystal grain).
2. The magnitude of the OP is rotationally invariant.

See this user forum post for details: https://groups.google.com/g/freud-users/c/tw_UAIeKD-Y/m/qlUTbZW9AgAJ

[vyasr]

Agree, we should probably add documentation indicating that the recommended usage is np.mean(np.abs(hex.particle_order)) or so for most standard use cases. I've never really seen much usage of the phase; it's also worth noting that strictly speaking the value of the phase is dependent on the choice of reference vector.

[bdice]

@glotzerlab/freud-developers In email correspondence, @Plastikschuessel noted that the literature seems to favor the 1/n normalization instead of 1/k, where n is the number of bonds for a particle and k is the order parameter symmetry. Given that the current implementation does not match the documentation and there is clear precedent in the literature, I think it would be reasonable to treat this as a bug and change the behavior to 1/n in the 2.x line of freud. This would also make Hexatic match the behavior of Steinhardt, which is its closest counterpart. I plan to file a small PR to correct this and will note the behavior change in the changelog.

[vyasr]

"Clear precedent" seems like a bit of a jump from "seems to favor"... could we document some of the relevant literature on this thread before making such a decision? I would also be curious if @joaander has any opinions since he's made extensive use of (and I believe he wrote) the feature.

[bdice]

Some literature / evidence for 1/n:

• @joaander’s 2017 paper, which I consider a definitive source for freud’s desired implementation: https://link.aps.org/doi/10.1103/PhysRevX.7.021001
• This uses k=n=6 so it’s not clear but I believe the intent is to use n, from the commentary: https://link.aps.org/doi/10.1103/PhysRevE.87.042134
• virtually all papers describing Steinhardt use 1/n (I am not aware of any counterexamples)

[vyasr]

Yes, I agree on the Steinhardt front. 1/n is the normalization that I'm most familiar with as well, I just want to make sure we can keep track of why we make these decisions in the event that we find issues with them later. But that seems fine to me then!

[bdice]

I have implemented this change in #657. I have a related question below:

The Hexatic order parameter is implemented in HexaticTranslational.cc alongside the Translational order parameter (which shares similar math -- both order parameters are implemented as subclasses of HexaticTranslational). The only use of k in the Translational order parameter is for normalization. If we change the normalization from 1/k to 1/n for Hexatic, should we also consider doing the same change to 1/n for the Translational order parameter? (The behavior in #657 is unchanged for Translational.) If so, the k parameter would not be needed and the API would need to change (or we'd need to document that k is an unused parameter and warn if it is provided).

[joaander]

Per https://aip.scitation.org/doi/10.1063/1.4774084, the closest approximation to the weighted morphometric neighborhood is to use n=k nearest neighbors. The only exception to that in https://link.aps.org/doi/10.1103/PhysRevX.7.021001 was triangles with n=3 nearest neighbors and k=6 order. I don't recommend using a fixed cutoff radius to compute hexatic order parameters, but if one were to do so, then 1/n is the the only reasonable choice for normalization - consider the case where there are two points at the center of each particle.

Is HexaticTranslational.cc the implementation of freud.order.Translational? What does that compute? (it isn't documented)

[bdice]

Is HexaticTranslational.cc the implementation of freud.order.Translational? What does that compute? (it isn't documented)

@joaander Yes, that is the implementation. Our C++ code lacks docstrings, and the Python docstring for the Translational order parameter is missing its formula definition. Good catch. I will add that in #657. If you happen to know of a paper that has used the Translational order parameter, it would be helpful to cite a specific reference for it. It's a "vestigial" feature whose origin is unknown, but we've kept it because it seemed moderately useful.

[joaander]

@bdice What is the formula that freud.order.Translational computes?

[bdice]

@joaander Here's what I put in #657.

freud/freud/order.pyx

Lines 397 to 407 in 72d1b9e

	The translational order parameter is used to measure order in the bonds
	of 2D systems. The translational order parameter for a particle :math:`i`
	and its :math:`n` neighbors :math:`j` is given by a sum over the
	neighbors, treating the 2D vectors between each pair of particles as a
	complex number with real part corresponding to the x-component of the
	vector and imaginary part corresponding to the y-component of the vector,
	divided by a normalization constant :math:`k`:
	:math:`\psi\left( i \right) = \frac{1}{k} \sum_j^n x_{ij} + y_{ij} i`
	The translational order parameter as written above is **complex-valued**.

[vyasr]

In the past when I've looked I've found various (pretty specific) definitions of a translational order parameter in the literature, but I don't recall finding anything that matched this definition. IIRC this is something either Michael or Wenbo wrote, and I have no idea if there is a canonical source for this.

[joaander]

I'm not aware of any works using that order parameter. How should it be used?

[vyasr]

I have absolutely no idea. Based on when the method was first created and who created it, my best guess is that it's something Wenbo and Michael developed at some stage of analyzing the data that led to Wenbo's 2D symmetry paper, but I don't recall it actually being used in that paper so I can't say for sure.

[bdice]

I see two use cases for Translational: coloring particles by the magnitude of the complex value, which gives a sense of the displacement, or using a cyclic color map applied to the complex arg (angle of the displacement vector). But likewise, I don't know of published works that made use of it.

[vyasr]

So option (1) is average bond length per particle, which can be trivially calculated from a neighborlist with distances, while option (2) gives you something like an average bond direction. I guess the former could be useful, not really sure about the latter, but I'm not really sure either of them merit a separate method. I think they can both be computed in <5 lines of Python and it wouldn't even be slow.

During earlier rounds of freud cleanup (around version 0.8 or so) I considered whether to remove this feature along with things like the kspace module, but I left it in place since I wasn't sure what it did and didn't want to remove features credited to particular developers without being sure how useful they might be. Maybe it's time to revisit that decision?

[bdice]

I support deprecation in 2.4 and removal in 3.0. Then we could simplify the joint HexaticTranslational implementation to be “only Hexatic.” Let’s use #659 to continue this conversation since it is unrelated to the original issue.