diff --git a/.cirrus.yml b/.cirrus.yml index 426c108..50f8c61 100644 --- a/.cirrus.yml +++ b/.cirrus.yml @@ -1,6 +1,6 @@ build_and_store_wheels: &BUILD_AND_STORE_WHEELS install_cibuildwheel_script: - - python -m pip install cibuildwheel==2.14.1 + - python -m pip install cibuildwheel==2.17.0 run_cibuildwheel_script: - cibuildwheel wheels_artifacts: @@ -24,11 +24,9 @@ linux_aarch64_task: macos_arm64_task: name: Build macOS arm64 wheels. macos_instance: - image: ghcr.io/cirruslabs/macos-monterey-xcode:latest + image: ghcr.io/cirruslabs/macos-sonoma-xcode env: - PATH: /opt/homebrew/opt/python@3.10/bin:$PATH - CIBW_ARCHS_MACOS: arm64 + PATH: /opt/homebrew/opt/python@3.10/libexec/bin:$PATH install_pre_requirements_script: - brew install python@3.10 - - ln -s python3 /opt/homebrew/opt/python@3.10/bin/python <<: *BUILD_AND_STORE_WHEELS diff --git a/.github/workflows/test-linux-build.yml b/.github/workflows/test-linux-build.yml index 2862dc8..b092259 100644 --- a/.github/workflows/test-linux-build.yml +++ b/.github/workflows/test-linux-build.yml @@ -32,7 +32,7 @@ jobs: uses: actions/cache/restore@v4 with: path: CASMcode_global/dist - key: ${{ runner.os }}-libcasm-global-v2-0-3 + key: ${{ runner.os }}-libcasm-global-v2-0-4 - name: Install CASM dependencies run: | diff --git a/.github/workflows/test-linux-cxx-only.yml b/.github/workflows/test-linux-cxx-only.yml index b0ad0e7..00e8e98 100644 --- a/.github/workflows/test-linux-cxx-only.yml +++ b/.github/workflows/test-linux-cxx-only.yml @@ -32,7 +32,7 @@ jobs: uses: actions/cache/restore@v4 with: path: CASMcode_global/dist - key: ${{ runner.os }}-libcasm-global-v2-0-3 + key: ${{ runner.os }}-libcasm-global-v2-0-4 - name: Install CASM dependencies run: | diff --git a/.github/workflows/test-linux-dependencies.yml b/.github/workflows/test-linux-dependencies.yml index fa33b9e..d8be1f5 100644 --- a/.github/workflows/test-linux-dependencies.yml +++ b/.github/workflows/test-linux-dependencies.yml @@ -24,7 +24,7 @@ jobs: uses: actions/cache/restore@v4 with: path: CASMcode_global/dist - key: ${{ runner.os }}-libcasm-global-v2-0-3 + key: ${{ runner.os }}-libcasm-global-v2-0-4 - name: checkout libcasm-global if: steps.cache-libcasm-global-restore.outputs.cache-hit != 'true' @@ -32,7 +32,7 @@ jobs: with: repository: prisms-center/CASMcode_global path: CASMcode_global - ref: v2.0.3 + ref: v2.0.4 - name: make global if: steps.cache-libcasm-global-restore.outputs.cache-hit != 'true' diff --git a/.github/workflows/test-linux.yml b/.github/workflows/test-linux.yml index 0c08077..3838d7b 100644 --- a/.github/workflows/test-linux.yml +++ b/.github/workflows/test-linux.yml @@ -32,7 +32,7 @@ jobs: uses: actions/cache/restore@v4 with: path: CASMcode_global/dist - key: ${{ runner.os }}-libcasm-global-v2-0-3 + key: ${{ runner.os }}-libcasm-global-v2-0-4 - name: Install CASM dependencies run: | diff --git a/CHANGELOG.md b/CHANGELOG.md index e31bc17..c5f576c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Fix CASM::xtal::make_primitive, which was not copying unique_names. This also fixes the output of libcasm.xtal.make_primitive which was losing the occ_dof list as a result. - Fix JSON output of xtal::BasicStructure site label +- Changed pseudoinverse calculation for basis changes to `completeOrthogonalDecomposition().pseudoInverse()` ### Changed @@ -28,6 +29,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Add to libcasm.xtal.Prim: method labels, constructor parameter `labels` - Add to libcasm.xtal.Lattice: methods reciprocal, volume, lengths_and_angles, from_lengths_and_angles - Added `unit_cell`, `diagonal_only`, and `fixed_shape` parameters to libcasm.xtal.enumerate_superlattices. +- Add to libcasm.xtal: combine_structures, filter_structure_by_atom_info, make_structure_atom_info, and make_structure_from_atom_info ## [2.0a8] - 2023-11-15 diff --git a/build_requirements.txt b/build_requirements.txt index 5ebd446..50998bb 100644 --- a/build_requirements.txt +++ b/build_requirements.txt @@ -4,4 +4,4 @@ scikit-build cmake>=3.20 ninja pybind11>=2.6 -libcasm-global>=2.0.2 +libcasm-global>=2.0.4 diff --git a/include/casm/crystallography/DoFSet.hh b/include/casm/crystallography/DoFSet.hh index 5f6ac57..3420b1c 100644 --- a/include/casm/crystallography/DoFSet.hh +++ b/include/casm/crystallography/DoFSet.hh @@ -41,11 +41,7 @@ class DoFSet { assert(m_basis.rows() == this->traits().dim()); // TODO: This makes sense I think? if (basis().cols() > 0 && basis().rows() > 0) { - m_inv_basis = basis() - .transpose() - .colPivHouseholderQr() - .solve(Eigen::MatrixXd::Identity(dim(), dim())) - .transpose(); + m_inv_basis = basis().completeOrthogonalDecomposition().pseudoInverse(); } } @@ -224,11 +220,6 @@ struct SiteDoFSetIsEquivalent_f : private DoFSetIsEquivalent_f { } }; -/// Create the symmtery representation for going from one basis to another -Eigen::MatrixXd dofset_transformation_matrix(const Eigen::MatrixXd &from_basis, - const Eigen::MatrixXd &to_basis, - double tol); - } // namespace xtal } // namespace CASM diff --git a/pyproject.toml b/pyproject.toml index 4536be2..515204e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,13 +6,13 @@ requires = [ "cmake>=3.20", "ninja", "pybind11>=2.6", - "libcasm-global>=2.0.2", + "libcasm-global>=2.0.4", ] build-backend = "setuptools.build_meta" [project] name = "libcasm-xtal" -version = "2.0a8" +version = "2.0a9" authors = [ { name="CASM developers", email="casm-developers@lists.engr.ucsb.edu" }, ] @@ -28,7 +28,7 @@ classifiers = [ "Topic :: Scientific/Engineering", ] dependencies = [ - "libcasm-global>=2.0.2", + "libcasm-global>=2.0.4", "numpy", ] @@ -51,13 +51,13 @@ select = ["E", "F", "I"] known-first-party = ["libcasm", "casm"] [tool.cibuildwheel] -# Build for python3.9, 3.10, 3.11 -build = "cp39-* cp310-* cp311-*" +# Build for python3.9, 3.10, 3.11, 3.12 +build = "cp39-* cp310-* cp311-* cp312-*" # Build for cpython only skip = "pp*" # Need libcasm-global at wheel repair stage -before-build = "pip install libcasm-global>=2.0.2" +before-build = "pip install libcasm-global>=2.0.4" # Testing test-requires = "pytest pytest-datadir" diff --git a/python/doc/conf.py b/python/doc/conf.py index fe33f3d..c65c86f 100644 --- a/python/doc/conf.py +++ b/python/doc/conf.py @@ -3,7 +3,7 @@ # -- package specific configuration -- project = "libcasm-xtal" version = "2.0" # The short X.Y version. -release = "2.0a8" # The full version, including alpha/beta/rc tags. +release = "2.0a9" # The full version, including alpha/beta/rc tags. project_desc = "CASM Crystallography" logo_text = "libcasm-xtal" github_url = "https://github.com/prisms-center/CASMcode_crystallography/" diff --git a/python/libcasm/xtal/__init__.py b/python/libcasm/xtal/__init__.py index 9c0c8c0..08d5df7 100644 --- a/python/libcasm/xtal/__init__.py +++ b/python/libcasm/xtal/__init__.py @@ -1,10 +1,14 @@ """CASM Crystallography""" from ._methods import ( StructureAtomInfo, + combine_structures, + filter_structure_by_atom_info, make_canonical, make_crystal_point_group, make_factor_group, make_primitive, + make_structure_atom_info, + make_structure_from_atom_info, make_within, sort_structure_by_atom_coordinate_cart, sort_structure_by_atom_coordinate_frac, diff --git a/python/libcasm/xtal/_methods.py b/python/libcasm/xtal/_methods.py index 03fbf60..6c67a72 100644 --- a/python/libcasm/xtal/_methods.py +++ b/python/libcasm/xtal/_methods.py @@ -1,7 +1,7 @@ import functools import math from collections import namedtuple -from typing import Any, Callable, Union +from typing import Any, Callable, Optional, Union import numpy as np @@ -211,7 +211,7 @@ def __lt__(self, other): "StructureAtomInfo", ["atom_type", "atom_coordinate_frac", "atom_coordinate_cart", "atom_properties"], ) -""" A namedtuple, used to hold atom info when sorting atoms in a +""" A namedtuple, used to hold atom info when sorting, filtering, etc. atoms in a :class:`~_xtal.Structure`. .. rubric:: Constructor @@ -228,10 +228,104 @@ def __lt__(self, other): :func:`~_xtal.Structure.atom_type.atom_coordinate_cart`. atom_properties: dict[str, numpy.ndarray[numpy.float64[m]]] The continuous properties associated with the atoms, if present, from - :func:`~_xtal.Structure.atom_type.atom_coordinate_frac`. + :func:`~_xtal.Structure.atom_type.atom_coordinate_frac`. All atoms + must have the same properties with values of the same dimension. """ +def make_structure_atom_info( + structure: _xtal.Structure, +) -> list[StructureAtomInfo]: + """Create a list of StructureAtomInfo from a Structure + + Parameters + ---------- + structure: _xtal.Structure + The structure to be sorted, filtered, etc. by atom info. + + Returns + ------- + structure_atom_info: list[StructureAtomInfo] + A list of StructureAtomInfo. + + """ + + atom_type = structure.atom_type() + atom_coordinate_frac = structure.atom_coordinate_frac() + atom_coordinate_cart = structure.atom_coordinate_cart() + atom_properties = structure.atom_properties() + + atoms = [] + import copy + + for i in range(len(atom_type)): + atoms.append( + StructureAtomInfo( + copy.copy(atom_type[i]), + atom_coordinate_frac[:, i].copy(), + atom_coordinate_cart[:, i].copy(), + {key: atom_properties[key][:, i].copy() for key in atom_properties}, + ) + ) + + return atoms + + +def make_structure_from_atom_info( + lattice: _xtal.Lattice, + atoms: list[StructureAtomInfo], + global_properties: dict[str, np.ndarray[np.float64]] = {}, +) -> _xtal.Structure: + """Create a Structure from a list of StructureAtomInfo + + Parameters + ---------- + lattice: _xtal.Lattice] + The lattice for the resulting structure. + atoms: list[StructureAtomInfo] + A list of StructureAtomInfo. The Cartesian coordinates are used when setting + atom coordinates in the resulting structure. + global_properties: dict[str, numpy.ndarray[numpy.float64[m, n]]] = {} + Continuous properties associated with entire crystal, if present. Keys must be + the name of a CASM-supported property type. Values are (m, 1) arrays with + dimensions matching the standard dimension of the property type. + + Returns + ------- + structure: _xtal.Structure + The resulting structure + """ + + n_atoms = len(atoms) + + atom_type = [atom.atom_type for atom in atoms] + atom_coordinate_cart = np.zeros((3, n_atoms)) + atom_properties = {} + + for i, atom in enumerate(atoms): + if i == 0: + for key, value in atom.atom_properties.items(): + dim = value.shape[0] + atom_properties[key] = np.zeros((dim, n_atoms)) + + atom_coordinate_cart[:, i] = atom.atom_coordinate_cart + for key, value in atom.atom_properties.items(): + atom_properties[key][:, i] = atom.atom_properties[key] + + atom_coordinate_frac = _xtal.cartesian_to_fractional( + lattice=lattice, + coordinate_cart=atom_coordinate_cart, + ) + + return _xtal.Structure( + lattice=lattice, + atom_type=atom_type, + atom_coordinate_frac=atom_coordinate_frac, + atom_properties=atom_properties, + global_properties=global_properties, + ) + + def sort_structure_by_atom_info( structure: _xtal.Structure, key: Callable[[StructureAtomInfo], Any], @@ -265,42 +359,16 @@ def sort_structure_by_atom_info( raise ValueError( "Error: only atomic structures may be sorted using sort_by_atom_info" ) - atom_type = structure.atom_type() - atom_coordinate_frac = structure.atom_coordinate_frac() - atom_coordinate_cart = structure.atom_coordinate_cart() - atom_properties = structure.atom_properties() - - atoms = [] - import copy - - for i in range(len(atom_type)): - atoms.append( - StructureAtomInfo( - copy.copy(atom_type[i]), - atom_coordinate_frac[:, i].copy(), - atom_coordinate_cart[:, i].copy(), - {key: atom_properties[key][:, i].copy() for key in atom_properties}, - ) - ) + atoms = make_structure_atom_info(structure) atoms.sort(key=key, reverse=reverse) - for i, atom in enumerate(atoms): - atom_type[i] = atom[0] - atom_coordinate_frac[:, i] = atom[1] - for key in atom_properties: - atom_properties[key][:, i] = atom[2][key] - - sorted_struc = _xtal.Structure( + return make_structure_from_atom_info( lattice=structure.lattice(), - atom_type=atom_type, - atom_coordinate_frac=atom_coordinate_frac, - atom_properties=atom_properties, + atoms=atoms, global_properties=structure.global_properties(), ) - return sorted_struc - def sort_structure_by_atom_type( structure: _xtal.Structure, @@ -474,3 +542,78 @@ def substitute_structure_species( mol_properties=structure.mol_properties(), global_properties=structure.global_properties(), ) + + +def filter_structure_by_atom_info( + structure: _xtal.Structure, + filter: Callable[[StructureAtomInfo], Any], +) -> _xtal.Structure: + """Return a copy of a structure with atoms passing a filter function + + .. rubric:: Example usage + + .. code-block:: Python + + # Remove all atoms with z coordinate >= 2.0 + structure_without_Al = filter_structure_by_atom_info( + input_structure, + lambda atom_info: atom_info.atom_coordinate_cart[2] < 2.0, + ) + + Parameters + ---------- + structure: _xtal.Structure + The initial structure + filter: Callable[[StructureAtomInfo], Any] + A function of `StructureAtomInfo` which returns True for atoms that should be + kept and False for atoms that should be removed. + + Returns + ------- + structure_after_removals: _xtal.Structure + A copy of `structure` excluding atoms for which the `filter` function returns + False. + """ + return make_structure_from_atom_info( + lattice=structure.lattice(), + atoms=[x for x in make_structure_atom_info(structure) if filter(x)], + global_properties=structure.global_properties(), + ) + + +def combine_structures( + structures: list[_xtal.Structure], + lattice: Optional[_xtal.Lattice] = None, + global_properties: dict[str, np.ndarray[np.float64]] = {}, +) -> _xtal.Structure: + """Return a new structure which combines the species of all input structures + + Parameters + ---------- + structures: list[_xtal.Structure] + The structures to be combined. + lattice: Optional[_xtal.Lattice] = None + If not None, the lattice of the resulting structure. If None, the resulting + structure will use the lattice of the first structure. + global_properties: dict[str, numpy.ndarray[numpy.float64[m, n]]] = {} + Continuous properties associated with entire crystal, if present. Keys must be + the name of a CASM-supported property type. Values are (m, 1) arrays with + dimensions matching the standard dimension of the property type. + + Returns + ------- + combined_structure: _xtal.Structure + The resulting structure, which contains the species of all input structures + with Cartesian coordinates fixed to the same values as in the input structures. + All atom or molecule properties remain the same as in the input structures. + """ + atoms = [] + for structure in structures: + if lattice is None: + lattice = structure.lattice() + atoms += make_structure_atom_info(structure) + return make_structure_from_atom_info( + lattice=lattice, + atoms=atoms, + global_properties=global_properties, + ) diff --git a/python/pyproject.toml b/python/pyproject.toml index 2d250d8..faf36ec 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -4,6 +4,6 @@ requires = [ "setuptools", "wheel", "pybind11>=2.8.0", - "libcasm-global>=2.0.2", + "libcasm-global>=2.0.4", ] build-backend = "setuptools.build_meta" diff --git a/python/setup.py b/python/setup.py index 07d2aa0..e6be42e 100644 --- a/python/setup.py +++ b/python/setup.py @@ -1,6 +1,6 @@ import os -__version__ = "2.0a8" +__version__ = "2.0a9" # Available at setup time due to pyproject.toml from pybind11.setup_helpers import Pybind11Extension, build_ext @@ -71,7 +71,7 @@ name="libcasm-xtal", version=__version__, packages=["libcasm", "libcasm.xtal"], - install_requires=["pybind11", "libcasm-global>=2.0.2"], + install_requires=["pybind11", "libcasm-global>=2.0.4"], ext_modules=ext_modules, cmdclass={"build_ext": build_ext}, ) diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp index 4553a03..02d705e 100644 --- a/python/src/xtal.cpp +++ b/python/src/xtal.cpp @@ -45,11 +45,7 @@ using namespace CASM; namespace _xtal_impl { Eigen::MatrixXd pseudoinverse(Eigen::MatrixXd const &M) { - Index dim = M.rows(); - return M.transpose() - .colPivHouseholderQr() - .solve(Eigen::MatrixXd::Identity(dim, dim)) - .transpose(); + return M.completeOrthogonalDecomposition().pseudoInverse(); } } // namespace _xtal_impl @@ -913,6 +909,8 @@ std::vector make_equivalent_property_values( if (basis.cols() == 0) { basis = Eigen::MatrixXd::Identity(dim, dim); } + // B * x_after = M * B * x_before + // x_after = B_pinv * M * B * x_before Eigen::MatrixXd basis_pinv = _xtal_impl::pseudoinverse(basis); for (auto const &op : point_group) { Eigen::VectorXd x_standard = basis * x; @@ -1297,7 +1295,8 @@ PYBIND11_MODULE(_xtal, m) { of `p` is undefined. )pbdoc"); - m.def("make_canonical_lattice", &make_canonical_lattice, py::arg("lattice"), + m.def("make_canonical_lattice", &make_canonical_lattice, + py::arg("init_lattice"), R"pbdoc( Returns the canonical equivalent lattice diff --git a/python/tests/test_structure.py b/python/tests/test_structure.py index 5e4bfa0..fd9f5a4 100644 --- a/python/tests/test_structure.py +++ b/python/tests/test_structure.py @@ -713,3 +713,51 @@ def test_substitute_structure_species_2(example_structure_1): {"A": "C", "B": "D"}, ) assert s2.atom_type() == ["C", "C", "D", "D"] + + +def test_combine_structures(example_structure_1): + lattice = xtal.Lattice( + np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]).transpose() + ) + + atom_coordinate_frac = np.array( + [ + [0.0, 0.0, 0.0], + [0.0, 0.5, 0.5], + [0.5, 0.0, 0.5], + [0.5, 0.5, 0.0], + ] + ).transpose() + + atom_type = ["A", "B", "B", "B"] + + s1 = xtal.Structure( + lattice=lattice, + atom_coordinate_frac=atom_coordinate_frac, + atom_type=atom_type, + ) + + combined = xtal.combine_structures([s1]) + assert combined.atom_type() == atom_type + print(xtal.pretty_json(combined.to_dict(frac=False))) + + translation = xtal.SymOp( + matrix=np.eye(3), + translation=np.array([0.0, 0.0, 3.0]), + time_reversal=False, + ) + s2 = translation * s1 + print(xtal.pretty_json(s2.to_dict(frac=False))) + combined = xtal.combine_structures( + structures=[s1, s2], + lattice=xtal.Lattice( + np.array( + [ + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 2.0], + ] + ).transpose() + ), + ) + assert combined.atom_type() == atom_type * 2 diff --git a/setup.py b/setup.py index aebf442..e7d518d 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,7 @@ setup( name="libcasm-xtal", - version="2.0a8", + version="2.0a9", packages=["libcasm", "libcasm.xtal"], package_dir={"": "python"}, cmake_install_dir="python/libcasm", diff --git a/src/casm/crystallography/DoFSet.cc b/src/casm/crystallography/DoFSet.cc index c049df4..7610775 100644 --- a/src/casm/crystallography/DoFSet.cc +++ b/src/casm/crystallography/DoFSet.cc @@ -49,15 +49,6 @@ bool DoFSetIsEquivalent_f::operator()(const DoFSet &other_value) const { this->_basis_spans_same_space(other_value); } -Eigen::MatrixXd dofset_transformation_matrix(const Eigen::MatrixXd &from_basis, - const Eigen::MatrixXd &to_basis, - double tol) { - Eigen::MatrixXd U = from_basis.colPivHouseholderQr().solve(to_basis); - if (!(U.transpose() * U).eval().isIdentity(tol)) { - throw std::runtime_error("Cannot find orthogonal symmetry representation!"); - } - return U; -} } // namespace xtal } // namespace CASM diff --git a/src/casm/crystallography/StrainConverter.cc b/src/casm/crystallography/StrainConverter.cc index 8ed3c54..7a6de03 100644 --- a/src/casm/crystallography/StrainConverter.cc +++ b/src/casm/crystallography/StrainConverter.cc @@ -14,11 +14,7 @@ namespace xtal { StrainConverter::StrainConverter(std::string _metric, Eigen::MatrixXd const &_basis) : m_metric(_metric), m_basis(_basis) { - Index dim = m_basis.cols(); - m_basis_pinv = m_basis.transpose() - .colPivHouseholderQr() - .solve(Eigen::MatrixXd::Identity(dim, dim)) - .transpose(); + m_basis_pinv = m_basis.completeOrthogonalDecomposition().pseudoInverse(); } /// \brief Decompose deformation tensor, F, as Q*U