From 8c8df25c7faa8bb200d8b915cf8ea72dda64d1d1 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Sun, 28 Jul 2024 19:45:58 -0700 Subject: [PATCH 01/11] point group of a cluster --- python/libcasm/xtal/__init__.py | 3 + python/libcasm/xtal/cluster_point_group.py | 1072 ++++++++++++++++++++ 2 files changed, 1075 insertions(+) create mode 100644 python/libcasm/xtal/cluster_point_group.py diff --git a/python/libcasm/xtal/__init__.py b/python/libcasm/xtal/__init__.py index 08d5df7..c951334 100644 --- a/python/libcasm/xtal/__init__.py +++ b/python/libcasm/xtal/__init__.py @@ -1,4 +1,5 @@ """CASM Crystallography""" + from ._methods import ( StructureAtomInfo, combine_structures, @@ -59,3 +60,5 @@ min_periodic_displacement, pretty_json, ) + +from .cluster_point_group import cluster_point_group diff --git a/python/libcasm/xtal/cluster_point_group.py b/python/libcasm/xtal/cluster_point_group.py new file mode 100644 index 0000000..22f4fb3 --- /dev/null +++ b/python/libcasm/xtal/cluster_point_group.py @@ -0,0 +1,1072 @@ +"""Contains functions related to finding the point group of a cluster of atoms or a molecule""" + +import itertools +import numpy as np +import libcasm.xtal as xtal +from typing import List, Tuple, Union + + +def geometric_center_of_mass(molecule: xtal.Occupant) -> np.ndarray: + """Calculates the geometric center of mass of the molecule + + :math:`\\mathbf{N} \\rightarrow` Number of atoms in molecule \n + :math:`\\mathbf{r}_i \\rightarrow` coordinates of :math:`i_{\\mathbf{th}}` atom \n + :math:`\\mathbf{r}_{com} \\rightarrow` coordinates of geometric center of mass + + .. math:: + \\mathbf{r}_{com} = \\frac{\\sum^{\\mathbf{N}}_{i=1} \\mathbf{r}_i}{\\mathbf{N}} + + Parameters + ---------- + molecule : xtal.Occupant + Molecule for which you want the center of mass + + Returns + ------- + numpy.ndarray + Coordinates of geometric center of mass + + """ + cart_coords = [atom.coordinate() for atom in molecule.atoms()] + geometric_center_of_mass = np.zeros(cart_coords[0].shape) + + for coord in cart_coords: + geometric_center_of_mass += coord + + return geometric_center_of_mass / len(cart_coords) + + +def shift_origin_of_molecule_to_geometric_center_of_mass( + molecule: xtal.Occupant, +) -> xtal.Occupant: + """Shifts the molecule to geometric center of mass + + :math:`\\mathbf{N} \\rightarrow` Group of all atoms in a molecule \n + :math:`\\mathbf{r}_i \\rightarrow` coordinates of :math:`i_{\\mathbf{th}}` atom \n + :math:`\\mathbf{r}_{com}` coordinates of geometric center of mass + + .. math:: + \\Big\\{ \\mathbf{r}_i - \\mathbf{r}_{com} \\Big| \\; \\forall \\; i \\; \\in \\; \\mathbf{N} \\Big\\} + + Parameters + ---------- + molecule : xtal.Occupant + Molecule for which you want shift the origin to it's geometric center of mass + + Returns + ------- + xtal.Occupant + Molecule where the origin is it's geometric center of mass + + """ + center_of_mass = geometric_center_of_mass(molecule) + new_atoms = [ + xtal.AtomComponent( + name=atom.name(), + coordinate=(atom.coordinate() - center_of_mass), + properties=atom.properties(), + ) + for atom in molecule.atoms() + ] + + return xtal.Occupant(name=molecule.name(), atoms=new_atoms) + + +def projection_operator(array: np.ndarray) -> np.ndarray: + """Make a projection operator for a given numpy array + + .. math:: + \\mathbf{R}_{proj} = \\mathbf{R} \\mathbf{R}^{T} + + Parameters + ---------- + array : numpy.ndarray + Numpy vector for which you want to compute projection operator + + Returns + ------- + numpy.ndarray + Projection operator + + """ + + if len(array.shape) == 1: + array = np.transpose(array[np.newaxis]) + + return array @ np.transpose(array) + + +def is_molecule_planar(molecule: xtal.Occupant, tol: float = 1e-5) -> bool: + """Checks if the provided molecule is planar + + :math:`\\mathbf{N} \\rightarrow` Group of all atoms in a molecule \n + :math:`\\mathbf{r}_i \\rightarrow` coordinates of :math:`i_{\\mathbf{th}}` atom \n + :math:`\\mathbf{r}^i_{proj} \\rightarrow` Projection operator of :math:`i_{\\mathbf{th}}` atom + (:math:`\\mathbf{r}_i \\mathbf{r}^{T}`)\n + + If only one of the eigen values of :math:`\\sum^{\\mathbf{N}}_{i=1} \\mathbf{r}^{i}_{proj}` + is zero, then the molecule is planar + + Parameters + ---------- + molecule : .core.Structure.Molecule + Molecule for which you want to compute it's planarity + tol : float, optional + Tolerance used for comparisions. + + Returns + ------- + bool + Whether the molecule is planar or not + """ + molecule_shifted_to_geometric_center_of_mass = ( + shift_origin_of_molecule_to_geometric_center_of_mass(molecule) + ) + + summed_projection_opertors_of_coords = projection_operators_sum_of_list_of_arrays( + get_coord_list_from_molecule(molecule_shifted_to_geometric_center_of_mass) + ) + + eigen_vals, _ = np.linalg.eig(summed_projection_opertors_of_coords) + + if ( + len( + [ + eigen_val + for eigen_val in eigen_vals + if np.allclose(eigen_val, 0, tol, tol) + ] + ) + == 1 + ): + return True + + return False + + +def is_molecule_linear(molecule: xtal.Occupant, tol: float = 1e-5) -> bool: + """Checks if the provided molecule is linear + + :math:`\\mathbf{N} \\rightarrow` Group of all atoms in a molecule \n + :math:`\\mathbf{r}_i \\rightarrow` coordinates of :math:`i_{\\mathbf{th}}` atom \n + :math:`\\mathbf{r}^i_{proj} \\rightarrow` Projection operator of :math:`i_{\\mathbf{th}}` atom + (:math:`\\mathbf{r}_i \\mathbf{r}^{T}`)\n + + If only two of the eigen values of :math:`\\sum^{\\mathbf{N}}_{i=1} \\mathbf{r}^{i}_{proj}` + is zero, then the molecule is linear + + Parameters + ---------- + molecule : xtal.Occupant + Molecule for which you want to compute it's linearity + tol : float, optional + Tolerance used for comparisions + + Returns + ------- + bool + Whether the molecule is linear or not + """ + + molecule_shifted_to_geometric_center_of_mass = ( + shift_origin_of_molecule_to_geometric_center_of_mass(molecule) + ) + + summed_projection_opertors_of_coords = projection_operators_sum_of_list_of_arrays( + get_coord_list_from_molecule(molecule_shifted_to_geometric_center_of_mass) + ) + + eigen_vals, _ = np.linalg.eig(summed_projection_opertors_of_coords) + + if ( + len( + [ + eigen_val + for eigen_val in eigen_vals + if np.allclose(eigen_val, 0, tol, tol) + ] + ) + == 2 + ): + return True + + return False + + +def projection_operators_sum_of_list_of_arrays( + array_list: List[np.ndarray], +) -> np.ndarray: + """For a given list of arrays, compute the projection operator of each array + and compute it's sum + + :math:`\\mathbf{N} \\rightarrow` Total number of arrays in the given list \n + :math:`\\mathbf{R}_i \\rightarrow` :math:`i_{\\mathbf{th}}` array in the list \n + + This function returns + + .. math:: + \\sum^{\\mathbf{N}}_{i=1} \\mathbf{R} \\mathbf{R}^{T} + + Parameters + ---------- + array_list : List[numpy.ndarray] + List of arrays + + Returns + ------- + numpy.ndarray + Summation of projection operators of each array in the given list + + """ + summed_projection_operators_array = np.zeros( + (array_list[0].shape[0], array_list[0].shape[0]) + ) + + for array in array_list: + summed_projection_operators_array += projection_operator(array) + + return summed_projection_operators_array + + +def check_if_two_molecules_are_equivalent( + molecule1: List[Tuple[str, np.ndarray]], + molecule2: List[Tuple[str, np.ndarray]], + tol: float = 1e-5, +) -> bool: + """Checks if the given molecules are equivalent + The molecules are provided as a list of tuples where each tuple represents + an atom type and it's corresponding cartesian coordinates.\n + Compares the given molecules by checking if each atom of a molecule 1 is also + present in molecule 2. + + Parameters + ---------- + molecule1 : List[Tuple[str, numpy.ndarray]] + Molecule 1 + molecule2 : List[Tuple[str, numpy.ndarray]] + Molecule 2 + tol : float, optional + Tolerance used for comparisions + + Returns + ------- + bool + Whether the given molecules are equivalent or not + + """ + if len(molecule1) == len(molecule2): + for atom1 in molecule1: + if ( + len( + [ + True + for atom2 in molecule2 + if (atom1[0] == atom2[0]) + and np.allclose(atom1[1], atom2[1], tol, tol) + ] + ) + == 0 + ): + return False + + return True + + return False + + +def apply_transformation_matrix_to_molecule( + transformation_matrix: np.ndarray, molecule: List[Tuple[str, np.ndarray]] +) -> List[Tuple[str, np.ndarray]]: + """Applies the given transformation matrix to the molecule + + :math:`\\mathbf{N} \\rightarrow` Group of all atoms in the molecule \n + :math:`\\mathbf{T} \\rightarrow` Transformation matrix \n + :math:`\\mathbf{r}_i \\rightarrow` Cartesian coordinates of :math:`i_{\\mathbf{th}}` atom in the molecule \n + :math:`el_i \\rightarrow` Element name of :math:`i_{\\mathbf{th}}` atom in the molecule \n + + This function returns + + .. math:: + \\Big \\{ ( el_i, \\mathbf{T} \\, \\mathbf{r}_i ) \\Big| \\; \\forall \\; i \\; \\in \\; \\mathbf{N} \\Big \\} + + Parameters + ---------- + transformation_matrix : numpy.ndarray + The number of columns of the transformation matrix should match + the dimensions of the coordinate space + molecule : List[Tuple[str, numpy.ndarray]] + Molecule to which you want to apply the given transformation matrix + + Returns + ------- + List[Tuple[str, numpy.ndarray]] + Transformed molecule + + Raises + ------ + ValueError + If number of columns of the ``transformation_matrix`` don't match + with the dimensions of the coordinate space + + """ + if ( + transformation_matrix.shape[len(transformation_matrix.shape) - 1] + != molecule[0][1].shape[0] + ): + raise ValueError( + "Columns of the transformation matrix do not match with the dimensions of the coordinate space" + ) + + return [(element, transformation_matrix @ coord) for element, coord in molecule] + + +def check_if_transformation_matrix_is_symmetry_opertion( + transformation_matrix: np.ndarray, + molecule: List[Tuple[str, np.ndarray]], + tol: float = 1e-5, +) -> bool: + """Checks if the transformation matrix you obtained in the process + of making the point group is a symmetry opertaion of the provided molecule. + + :math:`\\mathbf{T} \\rightarrow` Transformation matrix\n + The given transformation matrix is a symmetry operation under the following conditions:\n + * :math:`\\mathbf{T}` should be an orthogonal matrix (:math:`\\mathbf{T} \\mathbf{T}^T = \\mathbf{I}`) + * Applying :math:`\\mathbf{T}` to the given molecule should result in an equivalent molecule + + Parameters + ---------- + transformation_matrix : numpy.ndarray + Transformation matrix + molecule : List[Tuple[str, numpy.ndarray]] + Molecule represented as list of tuple of element name and it's corresponding cartesian + coordinate + tol : float, optional + Tolerance used for matrix comparisions + + Returns + ------- + bool + Whether the given transformation matrix is a symmetry operation + + """ + + # check for orthogonility + if np.allclose( + projection_operator(transformation_matrix), + np.identity(transformation_matrix.shape[0]), + tol, + tol, + ): + + # check whether the transformation matrix transforms the molecule into another equivalent molecule + transformed_molecule = apply_transformation_matrix_to_molecule( + transformation_matrix, molecule + ) + + if check_if_two_molecules_are_equivalent(molecule, transformed_molecule, tol): + return True + + return False + + +def convert_list_of_atom_coords_to_numpy_array(coords: List[np.ndarray]) -> np.ndarray: + """Converts a given list of coordinates into a numpy array where each + column is atom coordinates + + Parameters + ---------- + coords : List[numpy.ndarray] + List of coordinates that you want to convert to a numpy array + + Returns + ------- + numpy.ndarray + A matrix where columns of the matrix are coordinates of atoms + + """ + return np.transpose(np.array(coords)) + + +def get_coords_of_a_given_type_in_molecule( + molecule: xtal.Occupant, atom_type: str +) -> List[np.ndarray]: + """Given a molecule, return a list of cartesian coordinates of molecule for a specified + atom type + + Parameters + ---------- + molecule : xtal.Occupant + Molecule + atom_type : str + Atom type for which you want to get the cartesian coordinates + + Returns + ------- + List[numpy.ndarray] + List of cartesian coordinates of atoms in a given molecule belonging to a specified type + + """ + return [atom.coordinate() for atom in molecule.atoms() if atom.name() == atom_type] + + +def get_type_of_atoms_in_a_molecule(molecule: xtal.Occupant) -> List[str]: + """Returns all the unique atom types in a given molecule + + Parameters + --------- + molecule: xtal.Occupant + Molecule + + Returns + ------- + List[str] + List of all atom types in a molecule (Does not contain repeats) + + """ + + return list(set([atom.name() for atom in molecule.atoms()])) + + +def get_list_of_atom_type_and_coord_from_molecule( + molecule: xtal.Occupant, +) -> List[Tuple[str, np.ndarray]]: + """Given a molecule, return a list of tuples where each tuple contains an + atom type and it's cartesian coordinates + + Parameters + ---------- + molecule : xtal.Occupant + Molecule + + Returns + ------- + List[Tuple[str, numpy.ndarray]] + List of tuples where each tuple contains atom type and it's corresponding cartesian coordinates + + """ + return [(atom.name(), atom.coordinate()) for atom in molecule.atoms()] + + +def get_coord_list_from_molecule(molecule: xtal.Occupant) -> List[np.ndarray]: + """Get list of cartesian coordinates of a given molecule + + Parameters + ---------- + molecule : xtal.Occupant + Molecule + + Returns + ------- + List[numpy.ndarray] + List of all cartesian coordinates of a given molecule + + """ + return [atom.coordinate() for atom in molecule.atoms()] + + +def check_if_the_molecule_has_inversion_symmetry( + molecule: xtal.Occupant, tol: float = 1e-5 +) -> bool: + """ + Checks if the given molecule has inversion symmetry + + Parameters + ---------- + molecule : Molecule + Molecule + tol : float, optional + Tolerance used for comparisions + + Returns + ------- + bool + True if molecule as inversion symmetry, else False + + """ + + molecule = get_list_of_atom_type_and_coord_from_molecule(molecule) + transformation_matrix = np.identity(molecule[0][1].shape[0]) * -1 + molecule_after_inveresion = apply_transformation_matrix_to_molecule( + transformation_matrix, molecule + ) + + return check_if_two_molecules_are_equivalent( + molecule, molecule_after_inveresion, tol + ) + + +def is_symop_unique(sym_op: np.ndarray, sym_ops: List[np.ndarray], tol: float): + """Checks if the given sym_op is contained in the list of sym_ops + + Parameters + ---------- + sym_op : numpy.ndarray + sym_op that needs to be checked whether it is present in sym_ops + sym_ops : List[numpy.ndarray] + List of sym_ops in which whether the given sym_op is present or not + tol : float + Tolerance used for matrix comparisions + + Returns + ------- + bool + True if sym_op exists in sym_ops, else False + + """ + for s in sym_ops: + if np.allclose(sym_op, s, tol, tol): + return False + return True + + +def get_initial_combination( + number_of_atoms_of_a_given_type: int, is_planar: bool +) -> List[int]: + """ + Not intended for public use. Only to make the code more readable. + Get an initial combination of integers which can be used a pivoting point to find various + transformation matrices. + + The initial combination of integers depends on different conditions:\n + * If the molecule is planar:\n + * initial combination is [0, 1] + * If the molecule is non-planar:\n + * initial combination is [0, 1] if the number of atoms of a given type is 2 + * initial combination is [0, 1, 2] if the number of atoms of a given type is greater than 2 + + Parameters + ---------- + number_of_atoms_of_a_given_type : int + Number of atoms in a molecule with a specified type + is_planar : bool + Planarity of molecule + + Returns + ------- + List[int] + Initial combination of integers + + Raises + ------ + RuntimeError + If number of atoms is less than zero it happens + + """ + if number_of_atoms_of_a_given_type == 2 and is_planar is False: + initial_combination = list(range(2)) + elif number_of_atoms_of_a_given_type > 2 and is_planar is False: + initial_combination = list(range(3)) + elif number_of_atoms_of_a_given_type >= 2 and is_planar is True: + initial_combination = list(range(2)) + else: + raise ValueError( + "You should not have come here. If you did, please know that I am not programmed" + "to handle cases where number of atoms is less than or equal to zero." + ) + return initial_combination + + +def get_all_permutations( + number_of_atoms_of_a_given_type: int, is_planar: bool +) -> List[int]: + """ + Not intended for public use. Only to make the code more readable. + Get all permutations of integers which can be used to find various transformation matrices. + + :math:`\\mathbf{L} \\rightarrow` List of all integers - [0, ``number_of_atoms_of_a_given_type``]\n + This function returns:\n + * If the molecule is planar:\n + * All permutations and combinations of integers of dimensions 2 ([:math:`x`, :math:`y`]) from :math:`\\mathbf{L}` + * If the molecule is non-planar:\n + * If the ``number_of_atoms_of_a_given_type`` is 2, all permutations and combinations of integers of dimensions 2 [:math:`x`, :math:`y`] from :math:`\\mathbf{L}`\n + * If the ``number_of_atoms_of_a_given_type`` is greater than 2, all permutations and combinations of integers of dimensions 3 [:math:`x`, :math:`y`, :math:`z`] from :math:`\\mathbf{L}` + + Parameters + ---------- + number_of_atoms_of_a_given_type : int + Number of atoms in a molecule with a specified type + is_planar : bool + Planarity of molecule + + Returns + ------- + List[int] + List of all permutations + + Raises + ------ + RuntimeError + If number of atoms is less than zero it happens + """ + if number_of_atoms_of_a_given_type == 2 and is_planar is False: + all_permutations = list( + itertools.permutations(range(number_of_atoms_of_a_given_type), 2) + ) + elif number_of_atoms_of_a_given_type > 2 and is_planar is False: + all_permutations = list( + itertools.permutations(range(number_of_atoms_of_a_given_type), 3) + ) + + elif number_of_atoms_of_a_given_type >= 2 and is_planar is True: + all_permutations = list( + itertools.permutations(range(number_of_atoms_of_a_given_type), 2) + ) + else: + raise RuntimeError( + "You should not have come here. If you did, please know that I am not programmed" + "to handle cases where number of atoms is less than or equal to zero." + ) + return all_permutations + + +def get_transformation_matrices( + initial_coords: np.ndarray, + all_permuted_coords: List[np.ndarray], + number_of_atoms_of_a_given_type: int, + is_planar: bool, +) -> List[np.ndarray]: + """ + Not intended for public use. Only to make the code more readable. + Returns transformation matrices which might be symmetry operations. All these matrices + are generated from permuting atoms of same type. + + :math:`\\mathbf{r}_{init} \\rightarrow` A numpy matrix where columns of the matrix are cartesian coordinates + of atoms corresponding to an initial combination of atoms\n + :math:`\\mathbf{r}_{i} \\rightarrow` A numpy matrix where columns of the matrix are cartesian coordinates of atoms + corresponding to :math:`i_{\\mathbf{th}}` permutation\n + :math:`n \\rightarrow` Group of all permutations\n + :math:`\\mathbf{L} \\rightarrow` List of all transformation matrices\n + + If number of atoms of a given type is 2 and the molecule is non-planar:\n + + .. math:: + L = \\Big\\{ ( \\mathbf{r}_{i} \\mathbf{r}_{init}^{T})(\\mathbf{r}_{init} \\mathbf{r}_{init} ^{T})^{-1} + \\Big| \\; \\forall \\; i \\in \\; n \\Big\\} + + If number of atoms of a given type is greater than 2 and the molecule is non-planar or if the number + of atoms of a given type is greater than equal to 2 and the molecule is planar:\n + + .. math:: + L = \\Big\\{ \\mathbf{r}_{i} \\mathbf{r}_{init}^{-1} + \\Big| \\; \\forall \\; i \\in \\; n \\Big\\} + + Parameters + ---------- + initial_coords : numpy.ndarray + A numpy matrix where columns of the matrix are cartesian coordinates of atoms corresponding to an initial combination + all_permuted_coords : List[numpy.ndarray] + List of all matrices where columns of the matrix are cartesian coordinates of atoms corresponding to all permutations + and combinations + number_of_atoms_of_a_given_type : int + Number of atoms in a molecule of a specified type + is_planar : bool + Planarity of the molecule + + Returns + ------- + List[numpy.ndarray] + List of transformation matrices + + """ + transformation_matrices = [] * len(all_permuted_coords) + # If number of atoms of a given type is 2, C is not a square matrix + # T = (C_prime * transpose(C)) * inverse(C * transpose(C)) = (L_prime) * L_matrix + if number_of_atoms_of_a_given_type == 2 and is_planar is False: + l_matrix = np.linalg.inv(initial_coords @ np.transpose(initial_coords)) + for c_prime in all_permuted_coords: + l_prime = c_prime @ np.transpose(initial_coords) + t = l_prime @ l_matrix + transformation_matrices.append(t) + + # If number of atoms of a given type is > 2, + # T = C_prime * inverse(C) + if (number_of_atoms_of_a_given_type > 2 and is_planar is False) or ( + number_of_atoms_of_a_given_type >= 2 and is_planar is True + ): + for c_prime in all_permuted_coords: + t = c_prime @ np.linalg.inv(initial_coords) + transformation_matrices.append(t) + + return transformation_matrices + + +def make_symops_from_same_atom_type_coords( + atom_coords_of_a_given_type: List[np.ndarray], + molecule: List[Tuple[str, np.ndarray]], + is_planar: bool, + tol: float, +) -> List[np.ndarray]: + """ + Not intended for public use. Only to make the code more readable. + Makes symmetry operations from coordinates of same type of atoms. + + If the molecule is planar ``atom_coords_of_a_given_type`` should be reduced + dimensions (:math:`2 \\times 1``) and it should be in a coordinate space + where the planar molecule is always perpendicular to the :math:`z` axis. \n + + If the molecule is planar ``molecule`` should be in a coordinate space where + it is always perpendicular to the :math:`z` axis. \n + + To acheive these two criterion please use ``transformation_matrix_to_eigen_vector_space`` and + ``transform_coords_into_evs_for_planar_molecule`` helper functions. + + Parameters + ---------- + atom_coords_of_a_given_type : List[numpy.ndarray] + molecule : List[Tuple[str, numpy.ndarray]] + Molecule. Only used for checking if the obtained transformation matrix is a symmetry operation + is_planar : bool + Planarity of molecule + tol : float + Tolerance used for comparisons + + Returns + ------- + List[numpy.ndarray] + List of symmetry operations + + """ + + number_of_atoms_of_a_given_type = len(atom_coords_of_a_given_type) + sym_ops = [] + + # If you have only one atom of a given type + # All you can deduce is that Identity will be a sym op + if number_of_atoms_of_a_given_type == 1: + identity_symop = np.identity(molecule[0][1].shape[0]) + if is_symop_unique(identity_symop, sym_ops, tol): + sym_ops.append(identity_symop) + return sym_ops + + # Make an initial (3 x n) coordinate matrix where columns of the matrix are coords of atoms + initial_combination = get_initial_combination( + number_of_atoms_of_a_given_type, is_planar + ) + initial_coords = convert_list_of_atom_coords_to_numpy_array( + atom_coords_of_a_given_type + )[:, initial_combination] + + # Now get all the possible permutations and combinations of (3xn) matrices + # where columns of the matrices are cartesian coordinates atoms of a given type + all_permutations = get_all_permutations(number_of_atoms_of_a_given_type, is_planar) + all_permuted_coords = [ + convert_list_of_atom_coords_to_numpy_array(atom_coords_of_a_given_type)[ + :, permutation + ] + for permutation in all_permutations + ] + + # Find a transformation_matrices which takes initial coordinate matrix to permuted coordinates + # matrix. These are possible symmetry operations. + possible_sym_ops = get_transformation_matrices( + initial_coords, + all_permuted_coords, + number_of_atoms_of_a_given_type, + is_planar, + ) + + # Now these transformation matrices can only be symmetry operations of the entire + # molecule if and only if the molecule remains unchanged when you apply this transformation matrix + # to the whole molecule + if is_planar is False: + for possible_op in possible_sym_ops: + + if check_if_transformation_matrix_is_symmetry_opertion( + possible_op, molecule, tol + ) and is_symop_unique(possible_op, sym_ops, tol): + sym_ops.append(possible_op) + + if is_planar is True: + # If the molecule is planar, you end up getting 2x2 transformation matrices, since you are + # working in a reduced dimensions. Now you convert these 2x2 matrices into 3x3 matrices by + # adding another eigen vector +z or -z to it, since you are working in a space where the + # molecule is always perpendicular to z axis (this you have to ensure while passing molecule + # argument to the function). You also have to ensure that atom_coords_of_a_given_type will + # be 2 dimensional if your molecule is known to be planar. + new_possible_sym_ops = [] + for possible_op in possible_sym_ops: + new_possible_sym_ops.append( + np.block([[possible_op[0, :], 0], [possible_op[1, :], 0], [0, 0, 1]]) + ) + new_possible_sym_ops.append( + np.block([[possible_op[0, :], 0], [possible_op[1, :], 0], [0, 0, -1]]) + ) + + for op in new_possible_sym_ops: + if check_if_transformation_matrix_is_symmetry_opertion( + op, molecule, tol + ) and is_symop_unique(op, sym_ops, tol): + sym_ops.append(op) + + return sym_ops + + +def transformation_matrix_to_eigen_vector_space( + molecule: xtal.Occupant, tol: float +) -> np.ndarray: + """ + Transformation matrix that transforms the given planar molecule into eigen vector space + such that it is perpendicular to :math:`z` axis. + + Parameters + ---------- + molecule: xtal.Occupant + Planar molecule + tol : float + Tolerance used for comparisons + + Returns + ------- + numpy.ndarray + A matrix which transforms the planar molecule into eigen vector space such that it is + perpendicular to :math:`z` axis. + + """ + summed_projection_operators = projection_operators_sum_of_list_of_arrays( + get_coord_list_from_molecule(molecule) + ) + + eigen_vals, eigen_vector_space = np.linalg.eig(summed_projection_operators) + + index_corresponding_to_zero_eigen_val = [ + i + for i, eigen_val in enumerate(eigen_vals) + if np.allclose(eigen_val, 0, tol, tol) + ][0] + + coords_with_atom_types = get_list_of_atom_type_and_coord_from_molecule(molecule) + + # And now rotate the molecule in eigen vector space, so that the z coordinates of atom + # is always zero + identity = np.identity(coords_with_atom_types[0][1].shape[0]) + permutation = ( + list(range(0, index_corresponding_to_zero_eigen_val)) + + list( + range( + index_corresponding_to_zero_eigen_val + 1, + coords_with_atom_types[0][1].shape[0], + ) + ) + + [index_corresponding_to_zero_eigen_val] + ) + rotation_matrix = np.transpose(identity[:, permutation]) + + return rotation_matrix @ np.transpose(eigen_vector_space) + + +def transform_coords_into_reduced_evs_for_planar_molecule( + atom_type: str, + atom_coords_of_a_given_type: List[np.ndarray], + transf_mat_to_evs: np.ndarray, +) -> List[np.ndarray]: + """ + Not intended for public use. Only to make the code more readable. + This particular function gives atom coordinates for a given atom type + in a reduced (:math:`2 \\times 1`) dimensions in a eigen vector space + where the planar molecule is perpendicular to :math:`z` axis. + + Should only be used when the molecule is planar + + Parameters + ---------- + atom_type : str + Atom type + atom_coords_of_a_given_type : List[numpy.ndarray] + Atom coordinates of a given atom type in cartesian space + transf_mat_to_evs : numpy.ndarray + Transformation matrix which transforms the molecule into + eigen vector space where the planar molecule is perpendicular to :math:`z` axis + + Returns + ------- + List[numpy.ndarray] + Reduced dimension coordinates of atoms in eigen vector space where + the planar molecule is perpendicular to :math:`z` axis. + + """ + atom_coords_of_a_given_type_with_atom_type = [ + (atom_type, coord) for coord in atom_coords_of_a_given_type + ] + + atom_coords_of_a_given_type_in_evs_with_z_nulled = ( + apply_transformation_matrix_to_molecule( + transf_mat_to_evs, + atom_coords_of_a_given_type_with_atom_type, + ) + ) + atom_coords_of_a_given_type_in_evs_with_z_axed = [ + coord[1][0:-1] for coord in atom_coords_of_a_given_type_in_evs_with_z_nulled + ] + + return atom_coords_of_a_given_type_in_evs_with_z_axed + + +def issue_warning_if_origin_is_not_center_of_geometry( + molecule: xtal.Occupant, tol: float +) -> None: + """If center of geometry of a molecule is not at origin, + issue a warning + + Parameters + ---------- + molecule : Molecule + Molecule that you are interested in + tol : float + Tolerance for comparisions + + Returns + ------- + None + + """ + if not np.allclose( + geometric_center_of_mass(molecule), + np.zeros(molecule.atoms()[0].coordinate().shape), + tol, + tol, + ): + import warnings + + warnings.warn( + "WARNING: Geometric center of mass of the molecule is not at origin. Shifting it to origin. If you are working on the molecule further please shift it's geometric center of mass to origin." + ) + return None + + +def cluster_point_group( + molecule: xtal.Occupant, tol: float = 1e-5 +) -> Union[str, List[np.ndarray]]: + """Computes the point group of a cluster of atoms or a molecule + + The algorithm to compute the point group of molecule is divided into various parts:\n + + 1. If the molecule is linear, there can be infinite symmetry operations possible in the direction of + long axis, hence this function only returns the name of the point group.\n + + 2. If the molecule is planar,\n + The molecule is first transformed into eigen vector space and also will be made sure it is always + perpendicular to the :math:`z` axis. All the coordinates of atoms are represented in a reduced dimensions (:math:`2 \\times 1`) + in this vector space. Now for a particular atom type, you pick two atom coordinates each of dimension :math:`2\\times1`, + and make an initial :math:`2 \\times 2` (let's call it :math:`\\mathbf{R}`) matrix where + the columns of this matrix are coordinates of the two atoms. Now you make multiple other similar :math:`2\\times2` matrices by + looking at all the possible permutations and combinations of coordinates of atoms for the atom type you picked initially ( + let's call one of these matrices :math:`\\mathbf{R}^\\prime`). Now you compute various possible transformation matrices + by :math:`\\mathbf{R}^\\prime \\mathbf{R}^{-1}`. Now the dimensions of these transformation matrices are extended + to :math:`3 \\times 3` by adding the third axis :math:`\\pm z`. You repeat this process for multiple atom types and collect + all the transformation matrices. Now you filter symmetry operations from it by checking if the transformation + matrices are orthogonal and that the molecule is invariant to it. If either of these two criterion are not met, + the transformation matrix is discarded. Now these symmetry operations are transformed back into cartesian vector space. + + 3. If the molecule is non-planar,\n + a similar procedure as outlined for planar molecule is followed except you work in the normal cartesian space and all + the coordinates of atoms are :math:`3\\times1` vectors and the transformation matrices are :math:`3 \\times 3` matrices. + + Parameters + ---------- + molecule : xtal.Occupant + cluster of atoms for which you want the point group + tol : float, optional + Tolerance used for comparisons. The results will be + very sensitive if the given molecule has lot of + numerical noise. + + Returns + ------- + Union[str, List[numpy.ndarray]] + List of all 3x3 orthornormal matrices that form + the point group of a given molecule or a string if the + molecule is linear + + """ + + # Shift the molecule to geometric center of mass + issue_warning_if_origin_is_not_center_of_geometry(molecule, tol) + molecule = shift_origin_of_molecule_to_geometric_center_of_mass(molecule) + + # check if the molecular is linear + if is_molecule_linear(molecule, tol): + + # If the molecule has inversion symmetry and linear it belongs to D-inf Eh group + if check_if_the_molecule_has_inversion_symmetry(molecule, tol): + return "D\u221Eh" + + # If the molecule has no inversion symmetry and linear it belongs to C-inf Ev group + return "C\u221Ev" + + point_group = [] + # Identity is always a sym op + point_group.append(np.identity(molecule.atoms()[0].coordinate().shape[0])) + + # Get different atom types in the molecule + atom_types = get_type_of_atoms_in_a_molecule(molecule) + # Get molecule as a list of tuple of atom types and their corresponding cartesian coordinates + atom_coords_along_with_atom_types = get_list_of_atom_type_and_coord_from_molecule( + molecule + ) + + # If the molecular is planar and not linear + if ( + is_molecule_planar(molecule, tol) is True + and is_molecule_linear(molecule, tol) is False + ): + + # Transformation matrix that makes z coordinate of planar molecule zero + # You work in the eigen vector space where the molecule is perpendicular to the z axis + transf_mat_to_evs = transformation_matrix_to_eigen_vector_space(molecule, tol) + atom_coords_with_atom_types_in_evs_with_z_nulled = ( + apply_transformation_matrix_to_molecule( + transf_mat_to_evs, + atom_coords_along_with_atom_types, + ) + ) + + for atom_type in atom_types: + atom_coords_of_a_given_type = get_coords_of_a_given_type_in_molecule( + molecule, atom_type + ) + + # Transform the coordinates of each atom into eigen vector space such that the molecule is perpendicular + # to z axis. The coordinates here are considered as 2x1 matrices + atom_coords_of_a_given_type_in_evs_with_z_axed = ( + transform_coords_into_reduced_evs_for_planar_molecule( + atom_type, atom_coords_of_a_given_type, transf_mat_to_evs + ) + ) + + # The sym ops obtained by this routine are in the eigen vector space where molecule is perpendicular to + # to z axis + sym_ops = make_symops_from_same_atom_type_coords( + atom_coords_of_a_given_type_in_evs_with_z_axed, + atom_coords_with_atom_types_in_evs_with_z_nulled, + True, + tol, + ) + + for op in sym_ops: + # convert back the symmetry operations into cartesian vector space + new_op_in_old_space = np.linalg.inv(transf_mat_to_evs) @ op + new_op_in_old_space = new_op_in_old_space @ transf_mat_to_evs + if is_symop_unique(new_op_in_old_space, point_group, tol): + point_group.append(new_op_in_old_space) + + return point_group + + # If non planar + # For each atom type + for atom_type in atom_types: + # Get coordinates of one atom type + atom_coords_of_a_given_type = get_coords_of_a_given_type_in_molecule( + molecule, atom_type + ) + + sym_ops = make_symops_from_same_atom_type_coords( + atom_coords_of_a_given_type, + atom_coords_along_with_atom_types, + False, + tol, + ) + for sym_op in sym_ops: + if is_symop_unique(sym_op, point_group, tol): + point_group.append(sym_op) + + return point_group From 1e5d7a5cc9bbb76b20b9a8999720e91983bf901c Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Sun, 28 Jul 2024 20:05:32 -0700 Subject: [PATCH 02/11] cluster point group: return list of xtal.Symops --- python/libcasm/xtal/cluster_point_group.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/python/libcasm/xtal/cluster_point_group.py b/python/libcasm/xtal/cluster_point_group.py index 22f4fb3..eb88a9c 100644 --- a/python/libcasm/xtal/cluster_point_group.py +++ b/python/libcasm/xtal/cluster_point_group.py @@ -988,10 +988,16 @@ def cluster_point_group( # If the molecule has inversion symmetry and linear it belongs to D-inf Eh group if check_if_the_molecule_has_inversion_symmetry(molecule, tol): - return "D\u221Eh" + raise NotImplementedError( + "Method not implemented to get point group of a linear molecule" + ) + # return "D\u221Eh" # If the molecule has no inversion symmetry and linear it belongs to C-inf Ev group - return "C\u221Ev" + raise NotImplementedError( + "Method not implemented to get point group of a linear molecule" + ) + # return "C\u221Ev" point_group = [] # Identity is always a sym op @@ -1069,4 +1075,9 @@ def cluster_point_group( if is_symop_unique(sym_op, point_group, tol): point_group.append(sym_op) - return point_group + cluster_point_group = [ + xtal.SymOp(matrix=cart_op, translation=np.zeros(3), time_reversal=False) + for cart_op in point_group + ] + + return cluster_point_group From f462ad826d7a461eb351568a9b63b03f1e2b354f Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Tue, 30 Jul 2024 11:59:43 -0700 Subject: [PATCH 03/11] default properties arg for xtal.AtomComponent --- python/src/xtal.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp index eef9630..8aea9a2 100644 --- a/python/src/xtal.cpp +++ b/python/src/xtal.cpp @@ -1592,7 +1592,9 @@ PYBIND11_MODULE(_xtal, m) { An atomic component of a molecular :class:`Occupant` )pbdoc") .def(py::init(&make_atom_position), py::arg("name"), - py::arg("coordinate"), py::arg("properties"), R"pbdoc( + py::arg("coordinate"), + py::arg("properties") = std::map{}, + R"pbdoc( .. rubric:: Constructor From 3c3420e6863028a435e6e6bf46078632b0452a73 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Tue, 30 Jul 2024 12:00:18 -0700 Subject: [PATCH 04/11] doc update --- python/src/xtal.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp index 8aea9a2..530dd47 100644 --- a/python/src/xtal.cpp +++ b/python/src/xtal.cpp @@ -1608,7 +1608,7 @@ PYBIND11_MODULE(_xtal, m) { Position of the atom, in Cartesian coordinates, relative to the basis site at which the occupant containing this atom is placed. - properties : dict[str, array_like] + properties : dict[str, array_like], default={} Fixed properties of the atom, such as magnetic sping or selective dynamics flags. Keys must be the name of a CASM-supported property type. Values are arrays with @@ -1703,8 +1703,7 @@ PYBIND11_MODULE(_xtal, m) { Parameters ---------- name : str - A \"chemical name\", which must be identical for occupants - to be found symmetrically equivalent. The names are case + A \"chemical name\", which must be identical for occupants to be found symmetrically equivalent. The names are case sensitive, and "Va" is reserved for vacancies. )pbdoc"); From 1489c17d4024743c5438a4585ba5ca6f884ed738 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Tue, 30 Jul 2024 12:01:36 -0700 Subject: [PATCH 05/11] doc update --- python/src/xtal.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp index 530dd47..f58923f 100644 --- a/python/src/xtal.cpp +++ b/python/src/xtal.cpp @@ -1703,7 +1703,8 @@ PYBIND11_MODULE(_xtal, m) { Parameters ---------- name : str - A \"chemical name\", which must be identical for occupants to be found symmetrically equivalent. The names are case + A \"chemical name\", which must be identical for occupants + to be found symmetrically equivalent. The names are case sensitive, and "Va" is reserved for vacancies. )pbdoc"); From beb125b5537edfe28ec56ca8e7561e9a41f5fcc9 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Tue, 30 Jul 2024 15:28:01 -0700 Subject: [PATCH 06/11] xtal.Occupant.from_xyz_string() with tests; also fix test_prim.py shared_datadir fixture --- python/src/xtal.cpp | 58 ++++++++++++++++++++++++++++++++++++++- python/tests/conftest.py | 33 +++++++++++++++++++++- python/tests/test_prim.py | 20 +++++++------- 3 files changed, 99 insertions(+), 12 deletions(-) diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp index f58923f..1160412 100644 --- a/python/src/xtal.cpp +++ b/python/src/xtal.cpp @@ -289,6 +289,58 @@ std::map get_molecule_properties( return result; } +xtal::Molecule make_molecule_from_xyz_string( + std::string xyz_string, bool divisible = false, + std::map properties = {}) { + std::istringstream xyz_string_stream(xyz_string); + std::vector xyz_lines; + + std::string first_line; + std::getline(xyz_string_stream, first_line); + // the first line in the xyz format will be number of atoms in the molecule + int number_of_atoms_in_molecule = std::stoi(first_line); + + std::string name_of_molecule; + // the second line in the xyz format will be the name of the molecule + std::getline(xyz_string_stream, name_of_molecule); + + // from 3rd line, the xyz fomat follows "atom_type cart_coord_x cart_coord_y + // cart_coord_z" notation + + std::vector atom_positions; + for (std::string line; std::getline(xyz_string_stream, line);) { + // create a string stream for each line + // and parse atom_type and cartesian coordinates + std::stringstream line_stream(line); + + std::string atom_type; + line_stream >> atom_type; + + std::string cart_x_str; + line_stream >> cart_x_str; + + double cart_coord_x = std::stod(cart_x_str); + + std::string cart_y_str; + line_stream >> cart_y_str; + double cart_coord_y = std::stod(cart_y_str); + + std::string cart_z_str; + line_stream >> cart_z_str; + double cart_coord_z = std::stod(cart_z_str); + + xtal::AtomPosition atom_position( + Eigen::Vector3d(cart_coord_x, cart_coord_y, cart_coord_z), atom_type); + + atom_positions.push_back(atom_position); + } + + xtal::Molecule molecule(name_of_molecule, atom_positions, divisible); + molecule.set_properties(make_species_properties(properties)); + + return molecule; +} + // Prim std::shared_ptr make_prim( @@ -1687,7 +1739,11 @@ PYBIND11_MODULE(_xtal, m) { .def("is_vacancy", &xtal::Molecule::is_vacancy, "True if occupant is a vacancy.") .def("is_atomic", &xtal::Molecule::is_atomic, - "True if occupant is a single isotropic atom or vacancy"); + "True if occupant is a single isotropic atom or vacancy") + .def_static( + "from_xyz_string", &make_molecule_from_xyz_string, + py::arg("xyz_string"), py::arg("divisible") = false, + py::arg("properties") = std::map{}); m.def("make_vacancy", &xtal::Molecule::make_vacancy, R"pbdoc( Construct a Occupant object representing a vacancy diff --git a/python/tests/conftest.py b/python/tests/conftest.py index e32ea13..c12fe84 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -1,10 +1,41 @@ -import numpy as np +import os +import sys import pytest +import shutil +import pathlib +import numpy as np import libcasm.xtal as xtal import libcasm.xtal.lattices as xtal_lattices +def _win32_longpath(path): + """ + Helper function to add the long path prefix for Windows, so that shutil.copytree + won't fail while working with paths with 255+ chars. + """ + if sys.platform == "win32": + # The use of os.path.normpath here is necessary since "the "\\?\" prefix + # to a path string tells the Windows APIs to disable all string parsing + # and to send the string that follows it straight to the file system". + # (See https://docs.microsoft.com/pt-br/windows/desktop/FileIO/naming-a-file) + return "\\\\?\\" + os.path.normpath(path) + else: + return path + + +@pytest.fixture(scope="session") +def session_shared_datadir(tmpdir_factory): + original_shared_path = pathlib.Path(os.path.realpath(__file__)).parent / "data" + session_temp_path = tmpdir_factory.mktemp("session_data") + shutil.copytree( + _win32_longpath(original_shared_path), + _win32_longpath(str(session_temp_path)), + dirs_exist_ok=True, + ) + return session_temp_path + + @pytest.fixture def tetragonal_lattice(): # Lattice vectors diff --git a/python/tests/test_prim.py b/python/tests/test_prim.py index 1254138..7477ddd 100644 --- a/python/tests/test_prim.py +++ b/python/tests/test_prim.py @@ -75,16 +75,16 @@ def check_lial_with_occ_dofs(prim_with_occ_dofs, occ_dofs): assert prim_with_occ_dofs.labels() == [-1] * 4 -def test_prim_from_poscar(shared_datadir): - poscar_path = os.path.join(shared_datadir, "lial.vasp") +def test_prim_from_poscar(session_shared_datadir): + poscar_path = os.path.join(session_shared_datadir, "lial.vasp") # with no occ dofs prim = xtal.Prim.from_poscar(poscar_path) check_lial(prim) -def test_prim_from_poscar_with_occ_dof(shared_datadir): - poscar_path = os.path.join(shared_datadir, "lial.vasp") +def test_prim_from_poscar_with_occ_dof(session_shared_datadir): + poscar_path = os.path.join(session_shared_datadir, "lial.vasp") # with occ dofs occ_dofs = [["Li", "Va"], ["Li"], ["Al", "Va"], ["Li", "Al"]] @@ -92,8 +92,8 @@ def test_prim_from_poscar_with_occ_dof(shared_datadir): check_lial_with_occ_dofs(prim, occ_dofs) -def test_prim_from_poscar_str(shared_datadir): - poscar_path = os.path.join(shared_datadir, "lial.vasp") +def test_prim_from_poscar_str(session_shared_datadir): + poscar_path = os.path.join(session_shared_datadir, "lial.vasp") with open(poscar_path, "r") as f: poscar_str = f.read() @@ -103,8 +103,8 @@ def test_prim_from_poscar_str(shared_datadir): check_lial(prim) -def test_prim_from_poscar_str_with_occ_dof(shared_datadir): - poscar_path = os.path.join(shared_datadir, "lial.vasp") +def test_prim_from_poscar_str_with_occ_dof(session_shared_datadir): + poscar_path = os.path.join(session_shared_datadir, "lial.vasp") with open(poscar_path, "r") as f: poscar_str = f.read() @@ -149,8 +149,8 @@ def test_prim_from_poscar_str_selectivedynamics(): assert prim.local_dof() == [[], [], [], []] -def test_prim_to_poscar_str(shared_datadir): - poscar_path = os.path.join(shared_datadir, "lial.vasp") +def test_prim_to_poscar_str(session_shared_datadir): + poscar_path = os.path.join(session_shared_datadir, "lial.vasp") prim = xtal.Prim.from_poscar(poscar_path) structure = xtal.Structure( lattice=prim.lattice(), From 0fd56db8d2a34dfaf67d3772efcc50d6e0de3ee6 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Tue, 30 Jul 2024 15:29:11 -0700 Subject: [PATCH 07/11] add test_occupants.py --- python/tests/test_occupants.py | 41 ++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 python/tests/test_occupants.py diff --git a/python/tests/test_occupants.py b/python/tests/test_occupants.py new file mode 100644 index 0000000..306ae72 --- /dev/null +++ b/python/tests/test_occupants.py @@ -0,0 +1,41 @@ +import os +import numpy as np +from libcasm.xtal import Occupant + + +def test_occupant_from_xyz(session_shared_datadir): + """Test making xtal.Occupant.from_xyz_string() + method + + Parameters + ---------- + shared_datadir : str + + Returns + ------- + None + + """ + + with open( + os.path.join(session_shared_datadir, "input_molecules", "water.xyz"), "r" + ) as f: + water_xyz_string = f.read() + + water = Occupant.from_xyz_string(water_xyz_string) + + assert water.name() == "Water" + + expected_atom_names = ["O", "H", "H"] + expected_atom_coords = [ + np.array([0.0, 0.0, 0.11779]), + np.array([0.0, 0.75545, -0.47116]), + np.array([0.0, -0.75545, -0.47116]), + ] + + for atom, expected_atom_name, expected_atom_coord in zip( + water.atoms(), expected_atom_names, expected_atom_coords + ): + + assert atom.name() == expected_atom_name + assert np.allclose(atom.coordinate(), expected_atom_coord) From 260cc003bbf8d170b75cdcf14a72ebf61a681a90 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Tue, 30 Jul 2024 15:33:03 -0700 Subject: [PATCH 08/11] add molecule input files --- python/tests/data/input_molecules/ammonia.xyz | 6 ++ python/tests/data/input_molecules/benzene.xyz | 14 +++++ python/tests/data/input_molecules/c60.xyz | 62 +++++++++++++++++++ python/tests/data/input_molecules/ethane.xyz | 10 +++ python/tests/data/input_molecules/square.xyz | 6 ++ python/tests/data/input_molecules/water.xyz | 5 ++ 6 files changed, 103 insertions(+) create mode 100644 python/tests/data/input_molecules/ammonia.xyz create mode 100644 python/tests/data/input_molecules/benzene.xyz create mode 100644 python/tests/data/input_molecules/c60.xyz create mode 100644 python/tests/data/input_molecules/ethane.xyz create mode 100644 python/tests/data/input_molecules/square.xyz create mode 100644 python/tests/data/input_molecules/water.xyz diff --git a/python/tests/data/input_molecules/ammonia.xyz b/python/tests/data/input_molecules/ammonia.xyz new file mode 100644 index 0000000..cd36ea7 --- /dev/null +++ b/python/tests/data/input_molecules/ammonia.xyz @@ -0,0 +1,6 @@ +4 +Ammonia +N 0.000000 0.000000 0.000000 +H 0.000000 0.000000 1.008000 +H 0.950353 0.000000 -0.336000 +H -0.475176 -0.823029 -0.336000 diff --git a/python/tests/data/input_molecules/benzene.xyz b/python/tests/data/input_molecules/benzene.xyz new file mode 100644 index 0000000..dbb8bfd --- /dev/null +++ b/python/tests/data/input_molecules/benzene.xyz @@ -0,0 +1,14 @@ +12 +Benzene +H 1.2194 -0.1652 2.1600 +H -1.2644 -0.0630 2.1393 +H -2.4836 0.1021 -0.0204 +H -1.2194 0.1652 -2.1599 +H 1.2641 0.0628 -2.1395 +H 2.4836 -0.1022 0.0205 +C 0.6825 -0.0924 1.2087 +C -0.7075 -0.0352 1.1973 +C -1.3898 0.0572 -0.0114 +C -0.6824 0.0925 -1.2088 +C 0.7075 0.0352 -1.1973 +C 1.3899 -0.0572 0.0114 diff --git a/python/tests/data/input_molecules/c60.xyz b/python/tests/data/input_molecules/c60.xyz new file mode 100644 index 0000000..c48cc05 --- /dev/null +++ b/python/tests/data/input_molecules/c60.xyz @@ -0,0 +1,62 @@ +60 +C60 + C 2.16650 0.59060 2.58740 + C 3.03780 0.17660 1.59180 + C 1.27860 -0.30980 3.16790 + C 3.01180 -1.14340 1.16540 + C 3.10340 -1.43350 -0.19300 + C 3.15030 1.21060 0.66820 + C 3.24280 0.91490 -0.68590 + C 3.21920 -0.40230 -1.12070 + C -0.43930 1.35270 3.12710 + C 0.43630 2.26180 2.55420 + C -0.02960 0.06330 3.43790 + C 1.74420 1.87900 2.28300 + C 2.35190 2.26760 1.09900 + C -0.26330 3.02680 1.63260 + C 0.33740 3.40540 0.43730 + C 1.65160 3.02780 0.17070 + C -2.09030 -0.82250 2.59550 + C -2.51110 0.46640 2.28540 + C -0.84490 -1.02520 3.17380 + C -1.68740 1.55330 2.55120 + C -1.58430 2.58580 1.63190 + C -3.23140 0.40610 1.10070 + C -3.12270 1.44100 0.17460 + C -2.29470 2.52910 0.43990 + C -0.49080 -2.91330 1.73650 + C -1.74300 -2.71240 1.16370 + C -0.03930 -2.06840 2.74530 + C -2.54860 -1.66500 1.59420 + C -3.26020 -0.91410 0.67010 + C -1.65430 -3.00610 -0.18970 + C -2.35420 -2.24390 -1.11700 + C -3.16430 -1.19490 -0.68780 + C 2.13640 -2.05530 1.73580 + C 1.68950 -2.90090 0.72930 + C 1.27850 -1.63660 2.74350 + C 0.36780 -3.33270 0.73020 + C -0.34400 -3.39040 -0.45940 + C 2.28890 -2.52500 -0.46400 + C 1.57900 -2.57180 -1.65800 + C 0.25600 -3.00540 -1.65310 + C -2.18280 -0.57830 -2.59790 + C -1.74800 -1.86940 -2.30830 + C -0.43850 -2.24690 -2.58450 + C -1.28150 0.31890 -3.16710 + C -2.15260 2.05450 -1.73780 + C -3.04850 1.15350 -1.18110 + C -3.06560 -0.16290 -1.61070 + C -1.26610 1.64070 -2.72710 + C 0.50390 2.93610 -1.74180 + C -0.37880 3.35610 -0.75130 + C -1.69430 2.91860 -0.74910 + C 0.05210 2.07300 -2.73550 + C 2.09760 0.83400 -2.60510 + C 2.55170 1.69230 -1.61070 + C 1.75890 2.74520 -1.18240 + C 0.84200 1.02060 -3.17860 + C 0.44610 -1.34950 -3.16610 + C 1.69830 -1.54850 -2.59080 + C 2.51840 -0.46230 -2.31710 + C 0.02180 -0.06450 -3.45850 diff --git a/python/tests/data/input_molecules/ethane.xyz b/python/tests/data/input_molecules/ethane.xyz new file mode 100644 index 0000000..964206c --- /dev/null +++ b/python/tests/data/input_molecules/ethane.xyz @@ -0,0 +1,10 @@ +8 +Ethane +H 1.1851 -0.0039 0.9875 +H 1.1669 0.8330 -0.5693 +H 1.1155 -0.9329 -0.5145 +H -1.1669 -0.8334 0.5687 +H -1.1157 0.9326 0.5151 +H -1.1850 0.0044 -0.9875 +C -0.7516 0.0225 0.0209 +C 0.7516 -0.0225 -0.0209 diff --git a/python/tests/data/input_molecules/square.xyz b/python/tests/data/input_molecules/square.xyz new file mode 100644 index 0000000..0c73ccb --- /dev/null +++ b/python/tests/data/input_molecules/square.xyz @@ -0,0 +1,6 @@ +4 +Square +O 1.00000 1.00000 0.00000 +O -1.00000 1.00000 0.00000 +O -1.00000 -1.00000 0.00000 +O 1.00000 -1.00000 0.00000 diff --git a/python/tests/data/input_molecules/water.xyz b/python/tests/data/input_molecules/water.xyz new file mode 100644 index 0000000..ff96fe6 --- /dev/null +++ b/python/tests/data/input_molecules/water.xyz @@ -0,0 +1,5 @@ +3 +Water +O 0.00000 0.00000 0.11779 +H 0.00000 0.75545 -0.47116 +H 0.00000 -0.75545 -0.47116 From de58327c03ee181f5fdb04d6b821e0f12aed5ab0 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Tue, 30 Jul 2024 15:45:55 -0700 Subject: [PATCH 09/11] add tests for cluster_point_group --- python/libcasm/xtal/cluster_point_group.py | 4 +- python/tests/data/input_molecules/c60.xyz | 62 ------------ python/tests/test_cluster_point_group.py | 111 +++++++++++++++++++++ 3 files changed, 112 insertions(+), 65 deletions(-) delete mode 100644 python/tests/data/input_molecules/c60.xyz create mode 100644 python/tests/test_cluster_point_group.py diff --git a/python/libcasm/xtal/cluster_point_group.py b/python/libcasm/xtal/cluster_point_group.py index eb88a9c..0003170 100644 --- a/python/libcasm/xtal/cluster_point_group.py +++ b/python/libcasm/xtal/cluster_point_group.py @@ -925,9 +925,7 @@ def issue_warning_if_origin_is_not_center_of_geometry( tol, tol, ): - import warnings - - warnings.warn( + print( "WARNING: Geometric center of mass of the molecule is not at origin. Shifting it to origin. If you are working on the molecule further please shift it's geometric center of mass to origin." ) return None diff --git a/python/tests/data/input_molecules/c60.xyz b/python/tests/data/input_molecules/c60.xyz deleted file mode 100644 index c48cc05..0000000 --- a/python/tests/data/input_molecules/c60.xyz +++ /dev/null @@ -1,62 +0,0 @@ -60 -C60 - C 2.16650 0.59060 2.58740 - C 3.03780 0.17660 1.59180 - C 1.27860 -0.30980 3.16790 - C 3.01180 -1.14340 1.16540 - C 3.10340 -1.43350 -0.19300 - C 3.15030 1.21060 0.66820 - C 3.24280 0.91490 -0.68590 - C 3.21920 -0.40230 -1.12070 - C -0.43930 1.35270 3.12710 - C 0.43630 2.26180 2.55420 - C -0.02960 0.06330 3.43790 - C 1.74420 1.87900 2.28300 - C 2.35190 2.26760 1.09900 - C -0.26330 3.02680 1.63260 - C 0.33740 3.40540 0.43730 - C 1.65160 3.02780 0.17070 - C -2.09030 -0.82250 2.59550 - C -2.51110 0.46640 2.28540 - C -0.84490 -1.02520 3.17380 - C -1.68740 1.55330 2.55120 - C -1.58430 2.58580 1.63190 - C -3.23140 0.40610 1.10070 - C -3.12270 1.44100 0.17460 - C -2.29470 2.52910 0.43990 - C -0.49080 -2.91330 1.73650 - C -1.74300 -2.71240 1.16370 - C -0.03930 -2.06840 2.74530 - C -2.54860 -1.66500 1.59420 - C -3.26020 -0.91410 0.67010 - C -1.65430 -3.00610 -0.18970 - C -2.35420 -2.24390 -1.11700 - C -3.16430 -1.19490 -0.68780 - C 2.13640 -2.05530 1.73580 - C 1.68950 -2.90090 0.72930 - C 1.27850 -1.63660 2.74350 - C 0.36780 -3.33270 0.73020 - C -0.34400 -3.39040 -0.45940 - C 2.28890 -2.52500 -0.46400 - C 1.57900 -2.57180 -1.65800 - C 0.25600 -3.00540 -1.65310 - C -2.18280 -0.57830 -2.59790 - C -1.74800 -1.86940 -2.30830 - C -0.43850 -2.24690 -2.58450 - C -1.28150 0.31890 -3.16710 - C -2.15260 2.05450 -1.73780 - C -3.04850 1.15350 -1.18110 - C -3.06560 -0.16290 -1.61070 - C -1.26610 1.64070 -2.72710 - C 0.50390 2.93610 -1.74180 - C -0.37880 3.35610 -0.75130 - C -1.69430 2.91860 -0.74910 - C 0.05210 2.07300 -2.73550 - C 2.09760 0.83400 -2.60510 - C 2.55170 1.69230 -1.61070 - C 1.75890 2.74520 -1.18240 - C 0.84200 1.02060 -3.17860 - C 0.44610 -1.34950 -3.16610 - C 1.69830 -1.54850 -2.59080 - C 2.51840 -0.46230 -2.31710 - C 0.02180 -0.06450 -3.45850 diff --git a/python/tests/test_cluster_point_group.py b/python/tests/test_cluster_point_group.py new file mode 100644 index 0000000..29dfffb --- /dev/null +++ b/python/tests/test_cluster_point_group.py @@ -0,0 +1,111 @@ +import os +import pytest +from typing import Tuple +import libcasm.xtal as xtal + + +def supply_molecules_and_expected_point_group_operations( + session_shared_datadir: str, molecule_name: str +) -> Tuple[xtal.Occupant, str, float]: + """Supplies a list of some general molecules on which + you can test the point group code + + Parameters + ---------- + pytest_root_dir : str + molecule_name : str + + Returns + ------- + Tuple[xtal.Occupant, int, tol] + xtal.Occupant object, number of point group operations expected, + tolerance of comparisions + + + """ + molecules_dir = os.path.join(session_shared_datadir, "input_molecules") + + # For ammonia, you will have 6 point group operations + if molecule_name == "ammonia": + with open(os.path.join(molecules_dir, "ammonia.xyz"), "r") as f: + xyz_string = f.read() + return ( + xtal.Occupant.from_xyz_string(xyz_string), + 6, + 1e-3, + ) + + # For benzene, you will have 24 point group operations + elif molecule_name == "benzene": + with open(os.path.join(molecules_dir, "benzene.xyz"), "r") as f: + xyz_string = f.read() + return ( + xtal.Occupant.from_xyz_string(xyz_string), + 24, + 1e-3, + ) + + # For ethane, you will have 12 point group operations + elif molecule_name == "ethane": + with open(os.path.join(molecules_dir, "ethane.xyz"), "r") as f: + xyz_string = f.read() + return ( + xtal.Occupant.from_xyz_string(xyz_string), + 12, + 1e-2, + ) + + # For square, you will have 16 point group operations + elif molecule_name == "square": + with open(os.path.join(molecules_dir, "square.xyz"), "r") as f: + xyz_string = f.read() + return ( + xtal.Occupant.from_xyz_string(xyz_string), + 16, + 1e-5, + ) + + # For water, you will have 4 point group operations + elif molecule_name == "water": + with open(os.path.join(molecules_dir, "water.xyz"), "r") as f: + xyz_string = f.read() + return ( + xtal.Occupant.from_xyz_string(xyz_string), + 4, + 1e-5, + ) + + else: + raise RuntimeError("Unknown molecule type") + + +@pytest.mark.parametrize( + "molecule", + [ + "ammonia", + "benzene", + "ethane", + "square", + "water", + ], +) +def test_point_group_molecule(session_shared_datadir: str, molecule: str): + """Tests the point group of molecule for the given code + Currently tests against the hard coded number of point + group operations a molecule should have. + + Parameters + ---------- + pytest_root_dir : str + molecule: str + """ + ( + occupant, + expected_number_of_pg_operations, + tol, + ) = supply_molecules_and_expected_point_group_operations( + session_shared_datadir, molecule + ) + + point_group = xtal.cluster_point_group(occupant, tol) + assert len(point_group) == expected_number_of_pg_operations From 87bba8a989dce8193c0e84eea5c95bd0e9e6be38 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Tue, 30 Jul 2024 16:00:33 -0700 Subject: [PATCH 10/11] added pybind doc for xtal.Occupant.from_xyz_string() --- python/src/xtal.cpp | 42 +++++++++++++++++++++--- python/tests/test_cluster_point_group.py | 2 +- 2 files changed, 39 insertions(+), 5 deletions(-) diff --git a/python/src/xtal.cpp b/python/src/xtal.cpp index 1160412..8adf0e9 100644 --- a/python/src/xtal.cpp +++ b/python/src/xtal.cpp @@ -290,7 +290,7 @@ std::map get_molecule_properties( } xtal::Molecule make_molecule_from_xyz_string( - std::string xyz_string, bool divisible = false, + std::string xyz_string, bool is_divisible = false, std::map properties = {}) { std::istringstream xyz_string_stream(xyz_string); std::vector xyz_lines; @@ -335,7 +335,7 @@ xtal::Molecule make_molecule_from_xyz_string( atom_positions.push_back(atom_position); } - xtal::Molecule molecule(name_of_molecule, atom_positions, divisible); + xtal::Molecule molecule(name_of_molecule, atom_positions, is_divisible); molecule.set_properties(make_species_properties(properties)); return molecule; @@ -1742,8 +1742,42 @@ PYBIND11_MODULE(_xtal, m) { "True if occupant is a single isotropic atom or vacancy") .def_static( "from_xyz_string", &make_molecule_from_xyz_string, - py::arg("xyz_string"), py::arg("divisible") = false, - py::arg("properties") = std::map{}); + py::arg("xyz_string"), py::arg("is_divisible") = false, + py::arg("properties") = std::map{}, + R"pbdoc( + Make xtal.Occupant by reading the atom types and coordinates + of individual atoms represented in the xyz format + + Parameters + ---------- + xyz_string: str + String containing the atom types and coordinates + of the molecule in the xyz format + + xyz format for molecules uses the following notation: + ``` + + comment line (will be used to set the name of the molecule) + + + ... + ``` + is_divisible: bool + If True, indicates an Occupant that may split into components + during kinetic Monte Carlo calculations. + properties : dict[str, array_like], default={} + Fixed properties of the occupant, such as magnetic + spin or selective dynamics flags. Keys must be the name of a + CASM-supported property type. Values are arrays with + dimensions matching the standard dimension of the property + type. + + Returns + ------- + occupant: xtal.Occupant + )pbdoc" + + ); m.def("make_vacancy", &xtal::Molecule::make_vacancy, R"pbdoc( Construct a Occupant object representing a vacancy diff --git a/python/tests/test_cluster_point_group.py b/python/tests/test_cluster_point_group.py index 29dfffb..013f5fe 100644 --- a/python/tests/test_cluster_point_group.py +++ b/python/tests/test_cluster_point_group.py @@ -12,7 +12,7 @@ def supply_molecules_and_expected_point_group_operations( Parameters ---------- - pytest_root_dir : str + session_shared_datadir : str molecule_name : str Returns From 9c513f267d4adcbee3beb6bd7fc5c66966246449 Mon Sep 17 00:00:00 2001 From: seshasaibehara Date: Fri, 9 Aug 2024 14:02:04 -0700 Subject: [PATCH 11/11] revert unnecessary changes --- python/tests/conftest.py | 33 +--------------- python/tests/test_cluster_point_group.py | 10 ++--- python/tests/test_occupants.py | 6 +-- python/tests/test_prim.py | 50 +++++++++++++++++++----- 4 files changed, 47 insertions(+), 52 deletions(-) diff --git a/python/tests/conftest.py b/python/tests/conftest.py index c12fe84..e32ea13 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -1,41 +1,10 @@ -import os -import sys -import pytest -import shutil -import pathlib import numpy as np +import pytest import libcasm.xtal as xtal import libcasm.xtal.lattices as xtal_lattices -def _win32_longpath(path): - """ - Helper function to add the long path prefix for Windows, so that shutil.copytree - won't fail while working with paths with 255+ chars. - """ - if sys.platform == "win32": - # The use of os.path.normpath here is necessary since "the "\\?\" prefix - # to a path string tells the Windows APIs to disable all string parsing - # and to send the string that follows it straight to the file system". - # (See https://docs.microsoft.com/pt-br/windows/desktop/FileIO/naming-a-file) - return "\\\\?\\" + os.path.normpath(path) - else: - return path - - -@pytest.fixture(scope="session") -def session_shared_datadir(tmpdir_factory): - original_shared_path = pathlib.Path(os.path.realpath(__file__)).parent / "data" - session_temp_path = tmpdir_factory.mktemp("session_data") - shutil.copytree( - _win32_longpath(original_shared_path), - _win32_longpath(str(session_temp_path)), - dirs_exist_ok=True, - ) - return session_temp_path - - @pytest.fixture def tetragonal_lattice(): # Lattice vectors diff --git a/python/tests/test_cluster_point_group.py b/python/tests/test_cluster_point_group.py index 013f5fe..a8be628 100644 --- a/python/tests/test_cluster_point_group.py +++ b/python/tests/test_cluster_point_group.py @@ -5,7 +5,7 @@ def supply_molecules_and_expected_point_group_operations( - session_shared_datadir: str, molecule_name: str + shared_datadir: str, molecule_name: str ) -> Tuple[xtal.Occupant, str, float]: """Supplies a list of some general molecules on which you can test the point group code @@ -23,7 +23,7 @@ def supply_molecules_and_expected_point_group_operations( """ - molecules_dir = os.path.join(session_shared_datadir, "input_molecules") + molecules_dir = os.path.join(shared_datadir, "input_molecules") # For ammonia, you will have 6 point group operations if molecule_name == "ammonia": @@ -89,7 +89,7 @@ def supply_molecules_and_expected_point_group_operations( "water", ], ) -def test_point_group_molecule(session_shared_datadir: str, molecule: str): +def test_point_group_molecule(shared_datadir: str, molecule: str): """Tests the point group of molecule for the given code Currently tests against the hard coded number of point group operations a molecule should have. @@ -103,9 +103,7 @@ def test_point_group_molecule(session_shared_datadir: str, molecule: str): occupant, expected_number_of_pg_operations, tol, - ) = supply_molecules_and_expected_point_group_operations( - session_shared_datadir, molecule - ) + ) = supply_molecules_and_expected_point_group_operations(shared_datadir, molecule) point_group = xtal.cluster_point_group(occupant, tol) assert len(point_group) == expected_number_of_pg_operations diff --git a/python/tests/test_occupants.py b/python/tests/test_occupants.py index 306ae72..14492a8 100644 --- a/python/tests/test_occupants.py +++ b/python/tests/test_occupants.py @@ -3,7 +3,7 @@ from libcasm.xtal import Occupant -def test_occupant_from_xyz(session_shared_datadir): +def test_occupant_from_xyz(shared_datadir): """Test making xtal.Occupant.from_xyz_string() method @@ -17,9 +17,7 @@ def test_occupant_from_xyz(session_shared_datadir): """ - with open( - os.path.join(session_shared_datadir, "input_molecules", "water.xyz"), "r" - ) as f: + with open(os.path.join(shared_datadir, "input_molecules", "water.xyz"), "r") as f: water_xyz_string = f.read() water = Occupant.from_xyz_string(water_xyz_string) diff --git a/python/tests/test_prim.py b/python/tests/test_prim.py index 7477ddd..d8f353a 100644 --- a/python/tests/test_prim.py +++ b/python/tests/test_prim.py @@ -75,16 +75,16 @@ def check_lial_with_occ_dofs(prim_with_occ_dofs, occ_dofs): assert prim_with_occ_dofs.labels() == [-1] * 4 -def test_prim_from_poscar(session_shared_datadir): - poscar_path = os.path.join(session_shared_datadir, "lial.vasp") +def test_prim_from_poscar(shared_datadir): + poscar_path = os.path.join(shared_datadir, "lial.vasp") # with no occ dofs prim = xtal.Prim.from_poscar(poscar_path) check_lial(prim) -def test_prim_from_poscar_with_occ_dof(session_shared_datadir): - poscar_path = os.path.join(session_shared_datadir, "lial.vasp") +def test_prim_from_poscar_with_occ_dof(shared_datadir): + poscar_path = os.path.join(shared_datadir, "lial.vasp") # with occ dofs occ_dofs = [["Li", "Va"], ["Li"], ["Al", "Va"], ["Li", "Al"]] @@ -92,8 +92,8 @@ def test_prim_from_poscar_with_occ_dof(session_shared_datadir): check_lial_with_occ_dofs(prim, occ_dofs) -def test_prim_from_poscar_str(session_shared_datadir): - poscar_path = os.path.join(session_shared_datadir, "lial.vasp") +def test_prim_from_poscar_str(shared_datadir): + poscar_path = os.path.join(shared_datadir, "lial.vasp") with open(poscar_path, "r") as f: poscar_str = f.read() @@ -103,8 +103,8 @@ def test_prim_from_poscar_str(session_shared_datadir): check_lial(prim) -def test_prim_from_poscar_str_with_occ_dof(session_shared_datadir): - poscar_path = os.path.join(session_shared_datadir, "lial.vasp") +def test_prim_from_poscar_str_with_occ_dof(shared_datadir): + poscar_path = os.path.join(shared_datadir, "lial.vasp") with open(poscar_path, "r") as f: poscar_str = f.read() @@ -149,8 +149,8 @@ def test_prim_from_poscar_str_selectivedynamics(): assert prim.local_dof() == [[], [], [], []] -def test_prim_to_poscar_str(session_shared_datadir): - poscar_path = os.path.join(session_shared_datadir, "lial.vasp") +def test_prim_to_poscar_str(shared_datadir): + poscar_path = os.path.join(shared_datadir, "lial.vasp") prim = xtal.Prim.from_poscar(poscar_path) structure = xtal.Structure( lattice=prim.lattice(), @@ -233,6 +233,23 @@ def test_is_same_prim(simple_cubic_1d_disp_prim, simple_cubic_binary_prim): assert xtal._xtal._is_same_prim(second, first) is False +def test_copy(simple_cubic_binary_prim): + import copy + + prim = simple_cubic_binary_prim + prim1 = prim.copy() + assert isinstance(prim1, xtal.Prim) + assert prim1 is not prim + + prim2 = copy.copy(prim) + assert isinstance(prim2, xtal.Prim) + assert prim2 is not prim + + prim3 = copy.deepcopy(prim) + assert isinstance(prim3, xtal.Prim) + assert prim3 is not prim + + def test_to_dict(simple_cubic_binary_va_disp_Hstrain_prim): prim = simple_cubic_binary_va_disp_Hstrain_prim @@ -291,6 +308,19 @@ def test_from_dict(): assert prim_global_dof[0].dofname() == "Hstrain" +def test_repr(simple_cubic_binary_va_disp_Hstrain_prim): + import io + from contextlib import redirect_stdout + + prim = simple_cubic_binary_va_disp_Hstrain_prim + + f = io.StringIO() + with redirect_stdout(f): + print(prim) + out = f.getvalue() + assert "basis" in out + + def test_to_json(simple_cubic_binary_va_disp_Hstrain_prim): prim = simple_cubic_binary_va_disp_Hstrain_prim