Inconsistent reciprocalspace ASU between gemmi and cctbx

[JBGreisman]

Summary

As @kmdalton mentioned in #28, we have been implementing a few tests to make sure that our vectorized implementations of certain spacegroup-related operations in reciprocalspaceship are correct. In testing some code that maps Miller indices to the reciprocalspace ASU, we noticed some differences between our "reference implementations" in cctbx and gemmi.

We iterate over spacegroups using sgtbx.space_group_symbol_iterator(), which has 530 spacegroup settings, and use the extended Hermann-Mauguin for instantiating the gemmi.SpaceGroup. We are finding inconsistencies for 47/530 settings. This list of 47 spacegroups can be generated using the example below:

Example

reproduce_sg_inconsistency.py

from cctbx import sgtbx
from cctbx import miller
import gemmi
import numpy as np

def build_emptyMTZ(sg):
    """Build empty gemmi.MTZ"""
    m = gemmi.Mtz()
    m.spacegroup = sg
    m.add_dataset('crystal')
    m.add_column('H', 'H')
    m.add_column('K', 'H')
    m.add_column('L', 'H')
    m.add_column('M/ISYM', 'Y')
    return m
    
def gemmi_map2asu(H, symbol):
    """Map HKLs to reciprocal space ASU using gemmi"""
    xhm = symbol.universal_hermann_mauguin()
    sg = gemmi.SpaceGroup(xhm)
    m = build_emptyMTZ(sg)
    data = np.append(H, np.zeros((len(H), 1)), 1)
    m.set_data(data)
    m.switch_to_asu_hkl()
    h = m.column_with_label('H').array.astype(int)
    k = m.column_with_label('K').array.astype(int)
    l = m.column_with_label('L').array.astype(int)
    return np.vstack([h, k, l]).T
    
def cctbx_map2asu(H, symbol):
    """Map HKLs to reciprocal space ASU using cctbx"""
    sg  = sgtbx.space_group(symbol)
    asu = sgtbx.reciprocal_space_asu(sg.type())
    H_asu = np.vstack([miller.asym_index(sg, asu, i.tolist()).h() for i in H])
    return H_asu
    
hmin,hmax = -5, 5
H = np.mgrid[hmin:hmax+1:1,hmin:hmax+1:1,hmin:hmax+1:1].reshape((3, -1)).T

for symbol in sgtbx.space_group_symbol_iterator():
    xhm  = symbol.universal_hermann_mauguin()
    H_cctbx = cctbx_map2asu(H, symbol)
    H_gemmi = gemmi_map2asu(H, symbol)
    if not np.array_equal(H_cctbx, H_gemmi):
        print(xhm)

Output

Collapsed Output

I 1 2 1
I 1 1 2
I 2 1 1
P 1 n 1
P 1 1 n
P n 1 1
I 1 m 1
I 1 1 m
I m 1 1
A 1 n 1
I 1 a 1
I 1 c 1
B 1 1 n
I 1 1 b
A 1 1 n
I 1 1 a
C n 1 1
I c 1 1
B n 1 1
I b 1 1
I 1 2/m 1
I 1 1 2/m
I 2/m 1 1
P 1 2/n 1
P 1 1 2/n
P 2/n 1 1
P 1 21/n 1
P 1 1 21/n
P 21/n 1 1
A 1 2/n 1
I 1 2/a 1
I 1 2/c 1
B 1 1 2/n
I 1 1 2/b
A 1 1 2/n
I 1 1 2/a
C 2/n 1 1
I 2/c 1 1
B 2/n 1 1
I 2/b 1 1
R 3 :R
R -3 :R
R 3 2 :R
R 3 m :R
R 3 c :R
R -3 m :R
R -3 c :R

Single HKL example:

symbol = sgtbx.space_group_symbols("I 1 2 1") 
hkl = np.array([[-5, -5, 1]])
print(cctbx_map2asu(hkl, symbol)) # [[-5  5  1]]
print(gemmi_map2asu(hkl, symbol)) # [[ 5  5 -1]]

Issue

Is this expected behavior, or should gemmi and sgtbx both be using the same conventions for the ASU?

[wojdyr]

It's wrong in gemmi.
Thanks for the script that tests everything. I should be able to fix it today.

[wojdyr]

With the test program it was easy to fix it - by trials and errors :-)
Thanks again!

[JBGreisman]

Awesome -- I installed these changes locally, and I can confirm that all spacegroups are passing. Thanks for the quick attention to this!