Discrepancy between Phase.from_cif and Phase(structure)

Description of the change

Hello,

there appears to be a discrepancy in the atom xyz positions between creating a Phase using from_cif and assigning a structure directly (Phase(structure)).

This appears to be because lattice alignment is perfomed twice when using from_cif, once in the method itself and once in the Phase.structure setter.

I have removed the alignment in Phase.from_cif, as it will be called in the Phase constructor and added some checks to the test below, which were previously failing with the double lattice alignment.

Progress of the PR

[-] Docstrings for all functions
[x] Unit tests with pytest for all lines
[x] Clean code style by running black via pre-commit

Minimal example of the bug fix or new feature

from tempfile import NamedTemporaryFile

from orix.crystal_map import Phase
from diffpy.structure import loadStructure
import numpy as np

# from test
file_contents = """
#======================================================================

# CRYSTAL DATA

#----------------------------------------------------------------------

data_VESTA_phase_1

_chemical_name_common                  ''
_cell_length_a                         15.50000
_cell_length_b                         4.05000
_cell_length_c                         6.74000
_cell_angle_alpha                      90
_cell_angle_beta                       105.30000
_cell_angle_gamma                      90
_space_group_name_H-M_alt              'C 2/m'
_space_group_IT_number                 12

loop_
_space_group_symop_operation_xyz
   'x, y, z'
   '-x, -y, -z'
   '-x, y, -z'
   'x, -y, z'
   'x+1/2, y+1/2, z'
   '-x+1/2, -y+1/2, -z'
   '-x+1/2, y+1/2, -z'
   'x+1/2, -y+1/2, z'

loop_
   _atom_site_label
   _atom_site_occupancy
   _atom_site_fract_x
   _atom_site_fract_y
   _atom_site_fract_z
   _atom_site_adp_type
   _atom_site_B_iso_or_equiv
   _atom_site_type_symbol
   Mg(1)      1.0     0.000000      0.000000      0.000000     Biso  1.000000 Mg
   Mg(2)      1.0     0.347000      0.000000      0.089000     Biso  1.000000 Mg
   Mg(3)      1.0     0.423000      0.000000      0.652000     Biso  1.000000 Mg
   Si(1)      1.0     0.054000      0.000000      0.649000     Biso  1.000000 Si
   Si(2)      1.0     0.190000      0.000000      0.224000     Biso  1.000000 Si
   Al         1.0     0.211000      0.000000      0.626000     Biso  1.000000 Al"
"""

with NamedTemporaryFile(suffix=".cif", mode="w+", dir=".", delete=False) as f:
    f.write(file_contents)

structure = loadStructure(f.name)
phase = Phase(structure=structure)
phase_cif = Phase.from_cif(f.name)

# before fix
np.allclose(phase.structure.xyz, phase_cif.structure.xyz)
False

np.allclose(phase.structure.xyz_cartn, phase_cif.structure.xyz_cartn)
False

For reviewers

[x] The PR title is short, concise, and will make sense 1 year later.
[ ] New functions are imported in corresponding __init__.py.
[x] New features, API changes, and deprecations are mentioned in the unreleased section in CHANGELOG.rst.
[x] Contributor(s) are listed correctly in __credits__ in orix/__init__.py and in .zenodo.json.

Hi @dasilvaale,

Well spotted! And great job creating a fix right away.

As you show, removing the extra call to Lattice.setLatBase() fixes the issue. And your fix seems to be the right one to me.

To record your contributions, would you be so kind as to add your name (or GitHub handle) in between Viljar and Alexander in the credits list in orix/__init__.py and in .zenodo.json?

https://github.com/pyxem/orix/blob/856491e341efc695da7b483fd108491a425f884f/orix/__init__.py#L18-L19

https://github.com/pyxem/orix/blob/856491e341efc695da7b483fd108491a425f884f/.zenodo.json#L45-L53

Since this is a bug fix, I think we should release a 0.12.1 patch release right away. Could you change the base branch to main?

And, is it OK with you if I push to your branch to prepare it for release (after we merge)?

This appears to be because lattice alignment is perfomed twice when using from_cif, once in the method itself and once in the Phase.structure setter.

Setting the lattice base twice is actually not the issue. The issue in from_cif() is that we change the lattice before we extract the old cartesian atom coordinates xyz_cartn. These are calculated from the fractional atom coordinates xyz and the lattice base A (diffpy.structure code). As spotted and fixed by @viljarjf in #469, we have to update the xyz_cartn when we change A, since they are defined with respect to the lattice base.

Building on your example, we can compare the two routes in Phase(structure=structure) and Phase.from_cif(structure).

Steps in Phase.structure.setter:

from copy import deepcopy

structure2 = deepcopy(structure)

# 1
A_old = structure2.lattice.base
A_new = _new_structure_matrix_from_alignment(A_old, x="a", z="c*")
# 2
r_old = structure2.xyz_cartn
# 3
structure2.lattice.setLatBase(A_new)
# 4
structure2.xyz_cartn = r_old

assert not np.allclose(structure.xyz, structure2.xyz)
# True
assert np.allclose(structure.xyz_cartn, structure2.xyz_cartn)
# True

Steps in Phase.from_cif():

structure3 = deepcopy(structure)

# 1
A_old1 = structure3.lattice.base
A_new1 = _new_structure_matrix_from_alignment(A_old1, x="a", z="c*")
# 2
structure3.lattice.setLatBase(A_new1)
# 3
A_old2 = structure3.lattice.base
A_new2 = _new_structure_matrix_from_alignment(A_old2, x="a", z="c*")
# 4
r_old = structure3.xyz_cartn
# 5
structure3.lattice.setLatBase(A_new2)
# 6
structure3.xyz_cartn = r_old

assert not np.allclose(structure.xyz, structure3.xyz)
# False
assert np.allclose(structure.xyz_cartn, structure3.xyz_cartn)
# False

# Structure matrix unchanged after setting it the second time
assert np.allclose(A_new1, A_new2)
# True

Fixing steps in Phase.from_cif() by extracting the atom coordinates prior to the first change of lattice base:

structure4 = deepcopy(structure)

# 1
A_old1 = structure4.lattice.base
A_new1 = _new_structure_matrix_from_alignment(A_old1, x="a", z="c*")
# 2
r_old = structure4.xyz_cartn
structure4.lattice.setLatBase(A_new1)
# 3
A_old2 = structure4.lattice.base
A_new2 = _new_structure_matrix_from_alignment(A_old2, x="a", z="c*")
# 4
#r_old = structure4.xyz_cartn
# 5
structure4.lattice.setLatBase(A_new2)
# 6
structure4.xyz_cartn = r_old

assert not np.allclose(structure.xyz, structure4.xyz)
# True
assert np.allclose(structure.xyz_cartn, structure4.xyz_cartn)
# True

Hi @hakonanes thanks for looking into this and thanks for the explanations. I'm glad the fix makes sense and am happy to contribute.

I have updated the branch to main and added my affiliations to the files you mentioned.

And yes, please feel free to push to this branch.

@dasilvaale, thank you for this contribution! It has now been published in orix 0.12.1 on both PyPI and Anaconda :rocket:

pyxem / orix