Make changing spacegroups easy

[kmdalton]

Right now there is no obvious way to change a DataSet instance to a different spacegroup which possibly has a different basisop and/or reciprocal asu. For unmerged data, I think something along the lines of the following is probably sufficient:

def change_spacegroup(self, new_sg, inplace=False):
    if not inplace:
        ds = self.copy()
    else:
        ds = self

    if not isinstance(new_sg, gemmi.SpaceGroup):
        new_sg = gemmi.SpaceGroup(new_sg)

    ds.apply_symop(ds.spacegroup.basisop.inverse(), inplace=True)
    ds.apply_symop(new_sg.basisop, inplace=True)

    return ds

For merged data, I think in addition to the basisop application, you would want to expand the asu to P1 and select the new asu.

IMHO there is enough nuance to this problem that we should provide a well tested method for this.

[JBGreisman]

I agree that this makes a lot of sense to implement. There are a lot of gotchas to changing spacegroup, so it makes sense to include a built-in method that handles this correctly. Also, the current intuitive method of changing spacegroup by just calling the spacegroup setter-method will not end up doing anything more than storing a new gemmi.SpaceGroup object.

We have to make sure the convention being used is to first apply basisop then apply symops. Otherwise, I agree that your method descriptions seem appropriate.

Just a note about DataSet.hkl_to_observed():

• This method currently is called when loading unmerged MTZ files in order to map the Miller indices to their observed indices based on the ISYM value. However, this method does not apply any basis ops, so it doesn't necessarily fully map Miller indices to the "P1" indices, depending on the spacegroup. Is this an appropriate convention, or should we also apply spacegroup.basisop.inverse() in this case?

[JBGreisman]

As we discussed offline, another gotcha that is worth documenting here is that such a method may also have to update the DataSet.cell attribute to shuffle around unit cell parameters.

For example, this would be true if one wanted to change spacegroup from P 21 1 1 to the reference setting P 1 21 1

[JBGreisman]

I am just linking a related gemmi issue here in case the discussion becomes useful for implementing this: https://github.com/project-gemmi/gemmi/issues/87

[wojdyr]

FYI. I added partly-working function:

Mtz.reindex(op: gemmi.Op, verbose: bool = False)

This is equivalent to program gemmi reindex, see also project-gemmi/gemmi#87

I'll improve it when I find out what doesn't work – I haven't tested it much. You're of course welcome to test it. If you write me when it gives wrong results I should be able to fix it.

In particular, I expect that it doesn't work correctly for reindexing that changes the volume of unit cell – I didn't have such a test case.

[kmdalton]

We'll have to think about how to test that. Maybe we can use sgtbx to make reference data?

More to the point, I think the hard part about reindexing from a user perspective is knowing what the correct reindexing op is. Every time I need to work this out for my projects, I end up staring confusedly at the international tables. Shouldn't it be possible to automate? Is there a way to write a function get_reindexing_op(old_spacegroup, new_spacegroup)?

[wojdyr]

PHENIX, according to docs, has program phenix.reindex. I tried to use equivalent cctbx command, for example:

cctbx.python -m iotbx.command_line.reindex 5wkd_phases.mtz change_of_basis='-l,k,h+l'

but I got RuntimeError: iotbx Internal Error: .... I didn't investigate it and used CCP4 pointless instead. But I tested it only on one example (C2 -> I2).

I think the reindexing operator between different settings of the same space group would be the change-of-basis operator (possibly inverted, or a combination of two c-o-b operators). These operators are tabulated in Gemmi. But how to determine reindexing operator between different space groups? And how to check if reindexing is possible / makes sense? I don't know. Anyway, I'd worry about it later on, after reindexing with given operator is sorted.

[JBGreisman]

PHENIX also has phenix.reflection_file_converter which appears to have a relevant --change_to_space_group flag. It seems to require --expand_to_p1 to be used in conjunction in order to change a spacegroup to a lower symmetry:

phenix.reflection_file_converter HEWL_sg96.mtz --expand_to_p1 --change_to_space_group=4  --mtz=HEWL_sg4.mtz

I'm not quite sure if this is a useful way to generate test data for this sort of operation, but it may at least point us in relevant directions in cctbx.