Python SymmOp.from_rotation_matrix_and_translation_vector方法代码示例

# 需要导入模块: from pymatgen.core.operations import SymmOp [as 别名]
# 或者: from pymatgen.core.operations.SymmOp import from_rotation_matrix_and_translation_vector [as 别名]
 def from_spacegroup_number(sgnum):
     datadir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'sg_data')
     filename = str(sgnum).zfill(3) + "*"
     files = sorted(glob.glob(os.path.join(datadir, filename)))
     with open(files[0], "r") as fid:
         symmops = []
         rots = []
         lines = fid.readlines()
         sgname = lines[0].strip()
         for i in xrange(1, len(lines)):
             toks = re.split(",", lines[i].strip())
             if len(toks) == 3:
                 rot = np.zeros((3, 3))
                 trans = [0, 0, 0]
                 for j in xrange(3):
                     tok = toks[j]
                     m = re.search("([\+\-]*)([xyz])", tok)
                     if m:
                         factor = -1 if m.group(1) == "-" else 1
                         loc = ord(m.group(2)) - 120
                         rot[j, loc] = factor
                         tok = re.sub("([\+\-]*)([xyz])", "", tok)
                         if tok.strip() != '':
                             trans[j] = eval(tok)
                     rots.append(rot)
                 symmops.append(SymmOp.from_rotation_matrix_and_translation_vector(rot, trans))
         return Spacegroup(sgname, sgnum, symmops)

# 需要导入模块: from pymatgen.core.operations import SymmOp [as 别名]
# 或者: from pymatgen.core.operations.SymmOp import from_rotation_matrix_and_translation_vector [as 别名]
 def get_symmetry_operations(self, cartesian=False):
     """
     Return symmetry operations as a list of SymmOp objects.
     By default returns fractional coord symmops.
     But cartesian can be returned too.
     """
     (rotation, translation) = self.get_symmetry()
     symmops = []
     for rot, trans in zip(rotation, translation):
         if cartesian:
             rot = np.dot(self._structure.lattice.md2c, np.dot(rot, self._structure.lattice.mc2d))
             trans = np.dot(self._structure.lattice.md2c, trans)
         symmops.append(SymmOp.from_rotation_matrix_and_translation_vector(rot, trans))
     return symmops

# 需要导入模块: from pymatgen.core.operations import SymmOp [as 别名]
# 或者: from pymatgen.core.operations.SymmOp import from_rotation_matrix_and_translation_vector [as 别名]
    def _test_rot(self, rot, origin, fixed, to_fit, tol_atoms, tol_atoms_plus):
        found_map = False
        mapping_op = None
        for site in fixed:
            logger.debug("Trying candidate rotation : \n" + str(rot))
            if site.species_and_occu == origin.species_and_occu:
                shift = site.coords
                op = SymmOp.from_rotation_matrix_and_translation_vector(rot.rotation_matrix, shift)
                nstruct = apply_operation(to_fit, op)
                correspondance = OrderedDict()
                all_match = True
                biggest_dist = 0
                # check to see if transformed struct matches fixed structure
                for trans in nstruct:
                    cands = fixed.get_sites_in_sphere(trans.coords, tol_atoms_plus)
                    if len(cands) == 0:
                        logger.debug("No candidates found1")
                        all_match = False
                        break
                    cands = sorted(cands, key=lambda a: a[1])
                    (closest, closest_dist) = cands[0]
                    if closest_dist > tol_atoms or closest.species_and_occu != trans.species_and_occu:
                        logger.debug("Closest dist too large! closest dist = {}".format(closest_dist))
                        all_match = False
                        break
                    correspondance[trans] = closest
                    if closest_dist > biggest_dist:
                        biggest_dist = closest_dist

                if not all_match:
                    continue

                if not are_sites_unique(correspondance.values(), False):
                    all_match = False
                else:
                    for k, v in correspondance.items():
                        logger.debug(str(k) + " fits on " + str(v))

                # now check to see if the converse is true -- do all of the
                # sites of fixed match up with a site in toFit
                # this used to not be here. This fixes a bug.
                logger.debug("Checking inverse mapping")
                inv_correspondance = OrderedDict()

                # it used to be fixed.getNumSites() != nStruct.getNumSites()
                # so only when the number of sites are different but it's
                # actually better to allways check the reverse. This
                # elimininates weird situations where two atoms fit to one (reduced in the
                # unit cell)
                for fixed_site in fixed:
                    cands = nstruct.get_sites_in_sphere(fixed_site.coords, tol_atoms_plus)
                    if len(cands) == 0:
                        logger.debug("Rejected because inverse mapping does not fit - Step 1")
                        all_match = False
                        break

                    cands = sorted(cands, key=lambda a: a[1])
                    (closest, closest_dist) = cands[0]

                    if closest_dist > tol_atoms or closest.species_and_occu != fixed_site.species_and_occu:
                        all_match = False
                        logger.debug("Rejected because inverse mapping does not fit - Step 2")
                        break
                    inv_correspondance[fixed_site] = closest

                if all_match:
                    if not are_sites_unique(inv_correspondance.values(), False):
                        all_match = False
                        logger.debug("Rejected because two atoms fit to the same site for the inverse")

                if all_match:
                    self.inv_correspondance = inv_correspondance
                    logger.debug("Correspondance for the inverse")
                    for k, v in inv_correspondance.items():
                        logger.debug("{} fits on {}".format(k, v))

                # The smallest correspondance array shouldn't have any equivalent sites
                if fixed.num_sites != to_fit.num_sites:
                    logger.debug("Testing sites unique")
                    if not are_sites_unique(correspondance.values()):
                        all_match = False
                        logger.debug("Rejected because the smallest correspondance array has equivallent sites")
                        break

                if all_match:
                    found_map = True
                    mapping_op = op
                    self.correspondance = correspondance
                    break

        return (found_map, mapping_op)

# 需要导入模块: from pymatgen.core.operations import SymmOp [as 别名]
# 或者: from pymatgen.core.operations.SymmOp import from_rotation_matrix_and_translation_vector [as 别名]
    def _get_candidate_rotations(self, origin, fixed, to_fit):
        tol_shear = self._tolerance_cell_misfit
        fixed_basis = fixed.lattice.matrix.transpose()
        # need to generate candidate rotations ...
        lengths = fixed.lattice.abc

        shells = []
        for i in range(3):
            dr = lengths[i] * math.sqrt(tol_shear / 2)
            shell = to_fit.get_neighbors_in_shell(origin.coords, lengths[i], dr)
            logger.debug("shell {} radius={} dr={}".format(i, lengths[i], dr))
            shell = filter(lambda x: x[0].species_and_occu == origin.species_and_occu, shell)
            shell = sorted(shell, key=lambda x: x[1])
            shells.append([site for (site, dist) in shell])
            logger.debug("No. in shell = {}".format(len(shells[-1])))
        # now generate candidate rotations
        cand_rot = {}  # Dict of SymmOp : float
        a = len(shells[0])
        b = len(shells[1])
        c = len(shells[2])
        total_rots = a * b * c
        if total_rots < self._max_rotations:
            logger.debug("Total rots = {}. Using all rotations.".format(total_rots))
            test_rotations = itertools.product(*shells)
        else:
            logger.info(
                "Total rots = {m} exceed max_rotations = {n}. Using {n} randomly selected rotations.".format(
                    m=total_rots, n=self._max_rotations
                )
            )

            def random_rot():
                considered_rots = []
                while len(considered_rots) < self._max_rotations:
                    (x, y, z) = [random.randint(0, i - 1) for i in [a, b, c]]
                    if (x, y, z) not in considered_rots:
                        considered_rots.append((x, y, z))
                        yield (shells[0][x], shells[1][y], shells[2][z])

            test_rotations = random_rot()

        for pool in test_rotations:
            # now, can a unitary transformation bring the cell vectors together
            cell_v = np.array([nn.coords - origin.coords for nn in pool]).transpose()
            det = np.linalg.det(cell_v)
            if abs(det) < 0.001 or abs(abs(det) - fixed.volume) > 0.01:
                continue
            rot = np.dot(fixed_basis, np.linalg.inv(cell_v))
            r = SymmOp.from_rotation_matrix_and_translation_vector(rot, np.array([0, 0, 0]))

            if r not in cand_rot:
                transf = r.rotation_matrix
                transf = np.dot(transf.transpose(), transf)
                transf = np.eye(3) if almost_identity(transf) else transf
                pbis = sqrt_matrix(transf)
                shear_inv = shear_invariant(pbis)
                if shear_inv < tol_shear:
                    cand_rot[r] = shear_inv
                else:
                    logging.debug("Shear {} exceeds tol of {}".format(shear_inv, tol_shear))

        return cand_rot

# 需要导入模块: from pymatgen.core.operations import SymmOp [as 别名]
# 或者: from pymatgen.core.operations.SymmOp import from_rotation_matrix_and_translation_vector [as 别名]
    def fit(self, a, b):
        """
        Compares two structures and give the possible affine mapping that
        transforms one into another.
        """
        biggest_dist = 0
        logger.debug("Structure a")
        logger.debug(str(a))
        logger.debug("Structure b")
        logger.debug(str(b))

        # Check composition first.  If compositions are not the same, do not need to fit further.
        if a.composition.reduced_formula != b.composition.reduced_formula or (
            (a.num_sites != b.num_sites) and not self._supercells_allowed
        ):
            logger.debug("Compositions do not match")
            return None

        logger.debug("Compositions match")

        # Fitting is done by matching sites in one structure (to_fit)
        # to the other (fixed).
        # We set the structure with fewer sites as fixed,
        # and scale the structures to the same density
        (fixed, to_fit) = self._scale_structures(a, b)
        # Defines the atom misfit tolerance
        tol_atoms = self._tolerance_atomic_misfit * (3 * 0.7405 * fixed.volume / (4 * math.pi * fixed.num_sites)) ** (
            1 / 3
        )
        logger.debug("Atomic misfit tolerance = %.4f" % (tol_atoms))
        tol_atoms_plus = 1.1 * tol_atoms

        max_sites = float("inf")
        # determine which type of sites to use for the mapping
        for sp in to_fit.species_and_occu:
            sp_sites = [site for site in to_fit if site.species_and_occu == sp]
            if len(sp_sites) < max_sites:
                fit_sites = sp_sites
                max_sites = len(sp_sites)

        # Set the arbitrary origin
        origin = fit_sites[0]
        logger.debug("Origin = " + str(origin))
        # now that candidate rotations have been found, shift origin of to_fit
        # because the symmetry operations will be applied from this origin

        oshift = SymmOp.from_rotation_matrix_and_translation_vector(np.eye(3), -origin.coords)
        shifted_to_fit = apply_operation(to_fit, oshift)
        found_map = False

        # This is cheating, but let's try a simple rotation first.  In many situations,
        # E.g., when a structure has been topotatically delithiated or substituted, you actually
        # can get a very fast answer without having to try all rotations.
        simple_rot = self._get_rot_matrix(fixed, to_fit)
        if simple_rot is not None:
            rot = SymmOp.from_rotation_matrix_and_translation_vector(simple_rot, np.array([0, 0, 0]))
            (found_map, mapping_op) = self._test_rot(rot, origin, fixed, shifted_to_fit, tol_atoms, tol_atoms_plus)

        if not found_map:  # If simple rotation matching does not work, we have to search and try all rotations.
            logger.debug("Identity matching failed. Finding candidate rotations.")
            # Get candidate rotations
            cand_rot = self._get_candidate_rotations(origin, fixed, to_fit)

            logger.debug(" FOUND {} candidate rotations ".format(len(cand_rot)))
            if len(cand_rot) == 0:
                logger.debug("No candidate rotations found, returning null. ")
                return None

            # sort the operations, the first ones are the ones with small shear
            # this assures us that we find the smallest cell misfit fits
            sorted_cand_rot = sorted(cand_rot.keys(), key=lambda r: cand_rot[r])

            for rot in sorted_cand_rot:
                (found_map, mapping_op) = self._test_rot(rot, origin, fixed, shifted_to_fit, tol_atoms, tol_atoms_plus)

                if found_map:
                    break
        else:
            logger.debug("Identity matching found.")

        logger.debug("Done testing all candidate rotations")

        self._atomic_misfit = biggest_dist / ((3 * 0.7405 * a.volume / (4 * math.pi * a.num_sites)) ** (1 / 3))
        if mapping_op != None:
            rot = mapping_op.rotation_matrix  # maps toFit to fixed
            p = sqrt_matrix(np.dot(rot.transpose(), rot))
            scale_matrix = np.eye(3) * self.scale
            newrot = np.dot(scale_matrix, rot)
            # we need to now make sure fitterdata.MappingOp maps b -> a and not
            # the other way around
            mshift = np.dot(rot, oshift.translation_vector)
            finaltranslation = mapping_op.translation_vector + mshift[0]
            composite_op = SymmOp.from_rotation_matrix_and_translation_vector(newrot, finaltranslation)
            self._mapping_op = composite_op if self.fixed_is_a else composite_op.inverse
            # self._mapping_op = mapping_op
            self._cell_misfit = shear_invariant(p)