rs-station / reciprocalspaceship

Tools for exploring reciprocal space
https://rs-station.github.io/reciprocalspaceship/
MIT License
28 stars 12 forks source link

Adding a warning to apply_symop() #134

Closed DHekstra closed 2 years ago

DHekstra commented 2 years ago

When converting from a centered spacegroup like C2 to a non-centered one like P1, it is natural to apply a symop like gemmi.Op("1/2*x-1/2*y,1/2*x+1/2*y,z"). When not careful with the input Miller indices, this can result in fractional Miller indices. It may be good to add a check that new indices are integers.

JBGreisman commented 2 years ago

Regarding changing spacegroups, we have a long-standing issue (#5) to address this -- exactly because of all the gotchas that can come up with regard to centering. It hasn't been a high priority for now and will require some diligence with regard to testing to make sure we get it all right.

I think that sort of check for fractional Miller indices could be a good idea though because it would cause issues in other places in rs.

kmdalton commented 2 years ago

This maybe should even raise an error. I'm fine with whatever you all decide.

JBGreisman commented 2 years ago

Just to expand on this, the current behavior when this arises is to coerce the result to an integer value. This can result in something that isn't the applied symmetry operation. This occurs because of the implementation using np.floor_divide() in rs.utils.apply_to_hkl():

https://github.com/Hekstra-Lab/reciprocalspaceship/blob/ae0e2f538f0d287a7003ce056a7a2613586d4fbf/reciprocalspaceship/utils/symop.py#L4-L20

I can see two solutions:

  1. Implement in rs.utils.apply_to_hkl() -- we can change this function to use np.divide() and raise a warning or error if the resulting array has non-integer parts. The checking for non-integer parts could be handled with something like not np.any(np.mod(Hnew, 1)), however the error-checking does run the risk of slowing down functions that heavily rely on apply_to_hkl() such as rs.utils.hkl_to_asu() (we test whether this really is a bottleneck though).
  2. Implement error-checking in DataSet.apply_symop() -- this will require modifying rs.utils.apply_to_hkl() to use np.divide, and letting other functions handle the testing of this. My hesitation here is that it could require re-implementing the same checks in different places. However, by construction, hkl_to_asu() won't produce fractional Miller indices, so these checks may only be needed in the more general methods like apply_symop()
JBGreisman commented 2 years ago

Looking back at this, one can detect gemmi.Op instances that won't run into this problem:

import gemmi

# Can't produce non-integer HKLs if provided integer HKLs
op = gemmi.Op("x,x-y,2*x-z")
print(((np.array(op.rot)/op.DEN)%1==0).all())       # prints True

# Can cause problems given integer HKL input
op = gemmi.Op("1/2*x-1/2*y,1/2*x+1/2*y,z")
print(((np.array(op.rot)/op.DEN)%1==0).all())       # prints False

The test if ((np.array(op.rot)/op.DEN)%1==0).all() can be used to go between the two operating modes for this function.