Wrong Wyckoff position assignment

[lan496]

Initially reported in CompRhys/aviary#96

Consider the following C2/c (No. 15) structure

_cell_length_a   6.67047583
_cell_length_b   10.79514875
_cell_length_c   9.91142379
_cell_angle_alpha   90.00000000
_cell_angle_beta   80.18039461
_cell_angle_gamma   90.00000000
loop_
 _atom_site_type_symbol
 _atom_site_label
 _atom_site_symmetry_multiplicity
 _atom_site_fract_x
 _atom_site_fract_y
 _atom_site_fract_z
 _atom_site_occupancy
  Ca  Ca0  1  0.00000000  0.59536646  0.25000000  1
  Ca  Ca1  1  0.00000000  0.40463354  0.75000000  1
  Ca  Ca2  1  0.50000000  0.09536646  0.25000000  1
  Ca  Ca3  1  0.50000000  0.90463354  0.75000000  1
  Fe  Fe4  1  0.00000000  0.00000000  0.00000000  1
  Fe  Fe5  1  0.00000000  0.00000000  0.50000000  1
  Fe  Fe6  1  0.50000000  0.50000000  -0.00000000  1
  Fe  Fe7  1  0.50000000  0.50000000  0.50000000  1

Spglib assigns Wyckoff positions: Ca (4e) and Fe (4a).
Moyopy assigns Wyckoff positions: Ca (4e) and Fe (4d).

Because Euclidean Wyckoff sets for C2/c are (ab)(cd)(e)(f), these Wyckoff sequences are not equivalent up to Euclidean normalizer.
https://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-normsets?from=wycksets&gnum=15

[lan496]

spglib's std_positions

[[0.5        0.90463354 0.25      ]   # 4e
 [0.5        0.09536646 0.75      ]   # 4e
 [0.         0.         0.        ]                 # 4a
 [0.         0.         0.5       ]                # 4a
 [0.         0.40463354 0.25      ]    # 4e (y=0.40463354)
 [0.         0.59536646 0.75      ]     # 4e
 [0.5        0.5        0.        ]               # 4a
 [0.5        0.5        0.5       ]]            # 4a

moyopy's std_cell

[[ 0.    -0.845  0.25 ]     # 4e(y=-0.845)
 [ 0.5   -0.345  0.25 ]    # 4e
 [-0.5   -0.655 -0.25 ]   # 4e
 [ 0.    -0.155 -0.25 ]     # 4e

 [ 0.25  -0.25   0.   ]      # 4d
 [ 0.75   0.25   0.   ]       # 4d
 [-0.25  -0.25  -0.5  ]   # 4d
 [ 0.25   0.25  -0.5  ]]   # 4d

[lan496]

We need to standardize basis vectors of monoclinic lattice so that alpha=gamma=pi/2, pi/2 <= beta|<=2*pi/3

[lan496]

I guess this Wyckoff position slipping occurs when Euclidean normalizers are enlarged by special metrics.
For the C2/c case, it corresponds to the conditions when cos(beta) is close to integer multiples of c/a or a/c.

import numpy as np
from spglib import get_symmetry_dataset
from moyopy import Cell, MoyoDataset

a = 8.0
b = 8.0
c = 10.0
print(c / a, a / c)

y_4e = 0.6
positions = [
    # Fe (4a)
    [0.0, 0.0, 0.0],
    [0.0, 0.0, 0.5],
    [0.5, 0.5, 0.0],
    [0.5, 0.5, 0.5],
    # Ca (4e)
    [0, y_4e, 0.25],
    [0, -y_4e, 0.75],
    [0.5, 0.5 + y_4e, 0.25],
    [0.5, 0.5 - y_4e, 0.75],
]
numbers = [0, 0, 0, 0, 1, 1, 1, 1]
symprec = 1e-4

# https://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-normsets?from=wycksets&gnum=15
# Wyckoff sequence should be a-e or b-e
for beta in [np.deg2rad(deg) for deg in np.linspace(30, 150, num=128, endpoint=True)]:
    # C2/c (No. 15)
    lattice = [
        [a, 0.0, 0.0],
        [0.0, b, 0.0],
        [c * np.cos(beta), 0.0, c * np.sin(beta)],
    ]

    moyopy_cell = Cell(
        basis=lattice,
        positions=positions,
        numbers=numbers,
    )
    moyopy_dataset = MoyoDataset(moyopy_cell, symprec=symprec)
    print(f"{np.rad2deg(beta):.3f}, {np.cos(beta):.3f}, {sorted(list(set(moyopy_dataset.wyckoffs)))}")
    # print(f"Moyopy: {moyopy_dataset.number} {sorted(list(set(moyopy_dataset.wyckoffs)))}")