diff --git a/docs/source/api/diffpy.cmistructure.rst b/docs/source/api/diffpy.cmistructure.rst index f1eeec4..5332b6a 100644 --- a/docs/source/api/diffpy.cmistructure.rst +++ b/docs/source/api/diffpy.cmistructure.rst @@ -66,3 +66,11 @@ diffpy.cmistructure.sgconstraints module :members: :undoc-members: :show-inheritance: + +diffpy.cmistructure.converters module +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. automodule:: diffpy.cmistructure.converters + :members: + :undoc-members: + :show-inheritance: diff --git a/news/structure-to-crystal.rst b/news/structure-to-crystal.rst new file mode 100644 index 0000000..0079af1 --- /dev/null +++ b/news/structure-to-crystal.rst @@ -0,0 +1,23 @@ +**Added:** + +* ``create_crystal_from_structure`` in the new ``diffpy.cmistructure.converters`` module, which creates a pyobjcryst ``Crystal`` from a diffpy.structure ``Structure`` and a given space group. It reduces the atoms to the asymmetric unit and raises a ``ValueError`` if the lattice or the atoms don't fit the space group. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/src/diffpy/cmistructure/converters.py b/src/diffpy/cmistructure/converters.py new file mode 100644 index 0000000..3d4a9f1 --- /dev/null +++ b/src/diffpy/cmistructure/converters.py @@ -0,0 +1,163 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Members of the diffpy community. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Converters between structure representations.""" + +import numpy +from pyobjcryst.atom import Atom +from pyobjcryst.crystal import Crystal +from pyobjcryst.scatteringpower import ScatteringPowerAtom + +from diffpy.structure.symmetryutilities import epsilon, equal_positions + +__all__ = ["create_crystal_from_structure"] + +_ANISOTROPIC_ADP_NAMES = ("B11", "B22", "B33", "B12", "B13", "B23") + + +def create_crystal_from_structure(structure, space_group): + """Create a pyobjcryst Crystal from a diffpy.structure Structure. + + A Structure holds the full unit cell and does not record its space + group, so the space group must be given. The atoms are reduced to the + asymmetric unit of that space group, keeping the first atom of each set + of symmetry-equivalent atoms. The lattice, elements, positions, + occupancies and displacement factors are copied, with B = 8 pi^2 U. + Each atom is named by its label, or, if it has none, by its element and + a count, such as La1 or O2. + + The space group is read by ObjCryst, which takes Hermann-Mauguin + symbols with or without spaces, such as "P b n m" or "Pbnm", and space + group numbers, such as 62. A number gives the standard setting, so 62 + is "P n m a", not "P b n m". For a space group with two origin choices, + ObjCryst uses origin choice 2 unless another is given after a colon, + such as "F d -3 m:1". Give "P1" to keep every atom of the cell. + + Parameters + ---------- + structure : diffpy.structure.Structure + The structure to convert, holding the full unit cell. + space_group : str or int + The space group of the structure, as a Hermann-Mauguin symbol or a + space group number. + + Returns + ------- + pyobjcryst.crystal.Crystal + The crystal holding the asymmetric unit of the structure. + + Raises + ------ + ValueError + If ObjCryst does not recognize `space_group`, or if the lattice or + the atoms of `structure` do not fit it. + """ + lattice = structure.lattice + cell_parameters = [ + lattice.a, + lattice.b, + lattice.c, + lattice.alpha, + lattice.beta, + lattice.gamma, + ] + try: + crystal = Crystal( + lattice.a, + lattice.b, + lattice.c, + numpy.radians(lattice.alpha), + numpy.radians(lattice.beta), + numpy.radians(lattice.gamma), + str(space_group), + ) + except ValueError: + raise ValueError( + f"{space_group!r} is not a space group symbol or number. " + "Please give a Hermann-Mauguin symbol, such as 'P b n m', or a " + "number from 1 to 230." + ) from None + # ObjCryst silently changes a cell that doesn't fit the space group + crystal_cell_parameters = [ + crystal.a, + crystal.b, + crystal.c, + numpy.degrees(crystal.alpha), + numpy.degrees(crystal.beta), + numpy.degrees(crystal.gamma), + ] + if not numpy.allclose(crystal_cell_parameters, cell_parameters): + raise ValueError( + _get_not_fitting_message(crystal.GetSpaceGroup(), "lattice") + ) + element_counts = {} + for atom in _get_asymmetric_unit(structure, crystal.GetSpaceGroup()): + element_counts[atom.element] = element_counts.get(atom.element, 0) + 1 + name = atom.label or f"{atom.element}{element_counts[atom.element]}" + scattering_power = ScatteringPowerAtom(name, atom.element) + scattering_power.SetBiso(atom.Bisoequiv) + if atom.anisotropy: + for adp_name in _ANISOTROPIC_ADP_NAMES: + setattr(scattering_power, adp_name, getattr(atom, adp_name)) + x, y, z = atom.xyz + crystal.AddScatteringPower(scattering_power) + crystal.AddScatterer( + Atom(x, y, z, name, scattering_power, atom.occupancy) + ) + return crystal + + +def _get_asymmetric_unit(structure, space_group): + # Keep the first atom of each set of symmetry-equivalent atoms, and + # check that the structure holds every atom of each set + remaining_atoms = list(structure) + asymmetric_unit = [] + while remaining_atoms: + atom = remaining_atoms.pop(0) + asymmetric_unit.append(atom) + equivalent_atoms = [atom] + x, y, z = atom.xyz + positions = space_group.GetAllSymmetrics(x, y, z, False, False, True) + for position in positions: + matched_index = _find_atom_index( + equivalent_atoms, atom.element, position + ) + if matched_index is not None: + continue + index = _find_atom_index(remaining_atoms, atom.element, position) + if index is None: + raise ValueError( + _get_not_fitting_message(space_group, "atoms") + ) + equivalent_atoms.append(remaining_atoms.pop(index)) + return asymmetric_unit + + +def _find_atom_index(atoms, element, xyz): + for index, atom in enumerate(atoms): + if atom.element == element and equal_positions(atom.xyz, xyz, epsilon): + return index + return None + + +def _get_not_fitting_message(space_group, structure_part): + return ( + f"Space group {space_group.GetName()} does not fit the " + f"{structure_part} of the structure. Please give the full unit " + "cell of the structure and its space group, with the setting and " + "origin choice, or 'P1' to keep every atom." + ) + + +# End of file diff --git a/tests/conftest.py b/tests/conftest.py index f322641..6e8fdf8 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -12,7 +12,7 @@ from pyobjcryst.scatteringpower import ScatteringPowerAtom from diffpy.cmistructure.diffpyparset import DiffpyStructureParSet -from diffpy.structure import Structure +from diffpy.structure import Lattice, Structure C60_XYZ = """\ 3.451266498 0.685000000 0.000000000 @@ -136,6 +136,54 @@ def lamno3_p1_structure(datafile): return structure +@pytest.fixture +def lamno3_asymmetric_unit_structure(lamno3_p1_structure): + """Return the four atoms of the LaMnO3 asymmetric unit in Pbnm, + without the rest of the cell, as a diffpy Structure.""" + structure = lamno3_p1_structure[[0, 4, 8, 12]] + return structure + + +@pytest.fixture +def silicon_origin_1_structure(): + """Return the full silicon unit cell, with an atom at the origin, as + a diffpy Structure.""" + structure = Structure(lattice=Lattice(5.431, 5.431, 5.431, 90, 90, 90)) + positions = [ + [0.0, 0.0, 0.0], + [0.0, 0.5, 0.5], + [0.5, 0.0, 0.5], + [0.5, 0.5, 0.0], + [0.25, 0.25, 0.25], + [0.25, 0.75, 0.75], + [0.75, 0.25, 0.75], + [0.75, 0.75, 0.25], + ] + for xyz in positions: + structure.add_new_atom("Si", xyz=xyz) + return structure + + +@pytest.fixture +def silicon_origin_2_structure(): + """Return the full silicon unit cell, with the origin midway between + two bonded atoms, as a diffpy Structure.""" + structure = Structure(lattice=Lattice(5.431, 5.431, 5.431, 90, 90, 90)) + positions = [ + [0.125, 0.125, 0.125], + [0.125, 0.625, 0.625], + [0.625, 0.125, 0.625], + [0.625, 0.625, 0.125], + [0.875, 0.875, 0.875], + [0.875, 0.375, 0.375], + [0.375, 0.875, 0.375], + [0.375, 0.375, 0.875], + ] + for xyz in positions: + structure.add_new_atom("Si", xyz=xyz) + return structure + + @pytest.fixture def lamno3_generator(lamno3_p1_structure): parameter_set = DiffpyStructureParSet("LaMnO3", lamno3_p1_structure) diff --git a/tests/test_converters.py b/tests/test_converters.py new file mode 100644 index 0000000..7bec67c --- /dev/null +++ b/tests/test_converters.py @@ -0,0 +1,268 @@ +#!/usr/bin/env python +############################################################################## +# +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Members of the diffpy community. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst for license information. +# +############################################################################## +"""Tests for diffpy.cmistructure package.""" + +import re + +import numpy as np +import pytest + +from diffpy.cmistructure.converters import create_crystal_from_structure +from diffpy.structure import Atom, Lattice, Structure + + +@pytest.mark.parametrize( + "space_group", + [ + # Test converting the full LaMnO3 cell to its asymmetric unit + # C1: A user gives the space group with spaces, "P b n m". + # Expected: The crystal holds the four atoms of the asymmetric unit. + "P b n m", + # C2: A user gives the space group without spaces, "Pbnm". + # Expected: The crystal is the same as in C1. + "Pbnm", + ], +) +def test_create_crystal_from_structure(lamno3_p1_structure, space_group): + crystal = create_crystal_from_structure(lamno3_p1_structure, space_group) + # The lattice and space group are those of the structure + actual_lattice = [ + crystal.a, + crystal.b, + crystal.c, + np.degrees(crystal.alpha), + np.degrees(crystal.beta), + np.degrees(crystal.gamma), + ] + expected_lattice = [5.486341, 5.619215, 7.628206, 90.0, 90.0, 90.0] + assert np.allclose(actual_lattice, expected_lattice) + actual_space_group = crystal.GetSpaceGroup().GetName() + expected_space_group = "P b n m" + assert actual_space_group == expected_space_group + # The first atom of each set of symmetry-equivalent atoms is kept + atoms = [crystal.GetScatterer(i) for i in range(crystal.GetNbScatterer())] + actual_names = [atom.GetName() for atom in atoms] + expected_names = ["La1", "Mn1", "O1", "O2"] + assert actual_names == expected_names + actual_elements = [atom.GetScatteringPower().GetSymbol() for atom in atoms] + expected_elements = ["La", "Mn", "O", "O"] + assert actual_elements == expected_elements + actual_xyz = [[atom.X, atom.Y, atom.Z] for atom in atoms] + expected_xyz = [ + [0.99609631, 0.03214940, 0.25], + [0.0, 0.5, 0.0], + [0.05957463, 0.49616399, 0.25], + [0.72005206, 0.28938726, 0.03111255], + ] + assert np.allclose(actual_xyz, expected_xyz) + actual_occupancies = [atom.Occupancy for atom in atoms] + expected_occupancies = [1.0, 1.0, 1.0, 1.0] + assert np.allclose(actual_occupancies, expected_occupancies) + actual_biso = [atom.GetScatteringPower().GetBiso() for atom in atoms] + expected_uiso = np.array([0.00253993, 0.00065337, 0.00082010, 0.00512371]) + expected_biso = 8 * np.pi**2 * expected_uiso + assert np.allclose(actual_biso, expected_biso) + + +def test_create_crystal_from_structure_in_p1(): + # Test converting a structure to a crystal in P1 + # C1: A user keeps their whole cell in P1. The copper atom has a label, + # a partial occupancy and anisotropic displacements, and the oxygen + # atom has none of these. + # Expected: Both atoms are kept, the copper atom under its label and the + # oxygen atom as O1, each with its occupancy and displacement factors. + copper = Atom( + "Cu", + xyz=[0.1, 0.2, 0.3], + label="Cu_surface", + occupancy=0.8, + anisotropy=True, + U=[[0.004, 0.001, 0.0], [0.001, 0.006, 0.0], [0.0, 0.0, 0.008]], + ) + oxygen = Atom("O", xyz=[0.4, 0.5, 0.6], Uisoequiv=0.005) + lattice = Lattice(3.0, 4.0, 5.0, 90, 95, 90) + structure = Structure([copper, oxygen], lattice=lattice) + crystal = create_crystal_from_structure(structure, "P1") + actual_lattice = [ + crystal.a, + crystal.b, + crystal.c, + np.degrees(crystal.alpha), + np.degrees(crystal.beta), + np.degrees(crystal.gamma), + ] + expected_lattice = [3.0, 4.0, 5.0, 90.0, 95.0, 90.0] + assert np.allclose(actual_lattice, expected_lattice) + atoms = [crystal.GetScatterer(i) for i in range(crystal.GetNbScatterer())] + actual_names = [atom.GetName() for atom in atoms] + expected_names = ["Cu_surface", "O1"] + assert actual_names == expected_names + actual_xyz = [[atom.X, atom.Y, atom.Z] for atom in atoms] + expected_xyz = [[0.1, 0.2, 0.3], [0.4, 0.5, 0.6]] + assert np.allclose(actual_xyz, expected_xyz) + actual_occupancies = [atom.Occupancy for atom in atoms] + expected_occupancies = [0.8, 1.0] + assert np.allclose(actual_occupancies, expected_occupancies) + # The copper atom is anisotropic and the oxygen atom isotropic + copper_power = atoms[0].GetScatteringPower() + oxygen_power = atoms[1].GetScatteringPower() + actual_isotropic = [copper_power.IsIsotropic(), oxygen_power.IsIsotropic()] + expected_isotropic = [False, True] + assert actual_isotropic == expected_isotropic + actual_copper_bij = [ + copper_power.B11, + copper_power.B22, + copper_power.B33, + copper_power.B12, + copper_power.B13, + copper_power.B23, + ] + expected_copper_uij = np.array([0.004, 0.006, 0.008, 0.001, 0.0, 0.0]) + expected_copper_bij = 8 * np.pi**2 * expected_copper_uij + assert np.allclose(actual_copper_bij, expected_copper_bij) + actual_oxygen_biso = oxygen_power.GetBiso() + expected_oxygen_biso = 8 * np.pi**2 * 0.005 + assert np.isclose(actual_oxygen_biso, expected_oxygen_biso) + + +@pytest.mark.parametrize( + "structure_fixture_name, space_group, expected_space_group, expected_xyz", + [ + # Test converting silicon in F d -3 m, which has two origin choices + # C1: A user's silicon cell has an atom at the origin, and they ask + # for origin choice 1, "F d -3 m:1". + # Expected: The crystal holds one silicon atom, at the origin. + ( + "silicon_origin_1_structure", + "F d -3 m:1", + "F d -3 m:1", + [0.0, 0.0, 0.0], + ), + # C2: A user's silicon cell has its origin midway between two atoms, + # and they give "Fd-3m" without an origin choice. + # Expected: ObjCryst's default, origin choice 2, is used, and the + # crystal holds one silicon atom, at (1/8, 1/8, 1/8). + ( + "silicon_origin_2_structure", + "Fd-3m", + "F d -3 m:2", + [0.125, 0.125, 0.125], + ), + # C3: A user gives the space group number, 227, instead of the + # symbol. + # Expected: The number is read as "F d -3 m" with ObjCryst's default + # origin choice 2, as in C2. + ( + "silicon_origin_2_structure", + 227, + "F d -3 m:2", + [0.125, 0.125, 0.125], + ), + # C4: A user gives the space group number as text, "227". + # Expected: The crystal is the same as in C3. + ( + "silicon_origin_2_structure", + "227", + "F d -3 m:2", + [0.125, 0.125, 0.125], + ), + ], +) +def test_create_crystal_from_structure_with_origin_choice( + request, + structure_fixture_name, + space_group, + expected_space_group, + expected_xyz, +): + silicon_structure = request.getfixturevalue(structure_fixture_name) + crystal = create_crystal_from_structure(silicon_structure, space_group) + actual_space_group = crystal.GetSpaceGroup().GetName() + assert actual_space_group == expected_space_group + actual_atom_count = crystal.GetNbScatterer() + expected_atom_count = 1 + assert actual_atom_count == expected_atom_count + atom = crystal.GetScatterer(0) + actual_xyz = [atom.X, atom.Y, atom.Z] + assert np.allclose(actual_xyz, expected_xyz) + + +@pytest.mark.parametrize( + "structure_fixture_name, space_group, objcryst_space_group, " + "structure_part", + [ + # Test giving a space group that does not fit the structure + # C1: A user gives a cubic space group for the orthorhombic LaMnO3 + # cell. + # Expected: A ValueError says the space group doesn't fit the + # lattice. + ("lamno3_p1_structure", "Pm-3m", "P m -3 m", "lattice"), + # C2: A user refined LaMnO3 in Pbnm but gives Pnma, another setting + # of the same space group. + # Expected: A ValueError says the space group doesn't fit the atoms. + ("lamno3_p1_structure", "Pnma", "P n m a", "atoms"), + # C3: A user gives Fd-3m for a silicon cell with an atom at the + # origin, which needs origin choice 1, not ObjCryst's default of 2. + # Expected: A ValueError says the space group doesn't fit the atoms, + # naming origin choice 2. + ("silicon_origin_1_structure", "Fd-3m", "F d -3 m:2", "atoms"), + # C4: A user gives only the atoms of the LaMnO3 asymmetric unit, + # copied from a paper, instead of the full cell. + # Expected: A ValueError says the space group doesn't fit the atoms. + ("lamno3_asymmetric_unit_structure", "Pbnm", "P b n m", "atoms"), + ], +) +def test_create_crystal_from_structure_with_space_group_not_fitting( + request, + structure_fixture_name, + space_group, + objcryst_space_group, + structure_part, +): + structure = request.getfixturevalue(structure_fixture_name) + expected_msg = ( + f"Space group {objcryst_space_group} does not fit the " + f"{structure_part} of the structure. Please give the full unit " + "cell of the structure and its space group, with the setting and " + "origin choice, or 'P1' to keep every atom." + ) + with pytest.raises(ValueError, match=re.escape(expected_msg)): + create_crystal_from_structure(structure, space_group) + + +@pytest.mark.parametrize( + "space_group, space_group_in_message", + [ + # Test giving a space group that ObjCryst does not recognize + # C1: A user mistypes the space group symbol as "Pbnm2". + # Expected: A ValueError names the symbol and says how to give a + # space group. + ("Pbnm2", "'Pbnm2'"), + # C2: A user gives 231, which is not a space group number. + # Expected: A ValueError names the number and says how to give a + # space group. + (231, "231"), + ], +) +def test_create_crystal_from_structure_with_unknown_space_group( + lamno3_p1_structure, space_group, space_group_in_message +): + expected_msg = ( + f"{space_group_in_message} is not a space group symbol or number. " + "Please give a Hermann-Mauguin symbol, such as 'P b n m', or a " + "number from 1 to 230." + ) + with pytest.raises(ValueError, match=re.escape(expected_msg)): + create_crystal_from_structure(lamno3_p1_structure, space_group)