Source code for moldesign.geom.symmetry

# Copyright 2016 Autodesk Inc.
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.

import numpy as np

import moldesign as mdt
from moldesign import units as u
from moldesign.external import transformations as trns
from moldesign.interfaces import symmol_interface as smi

get_symmetry = smi.run_symmol

[docs]class SymmetryElement(object): """ Represents an basic, origin-centered symmetry operation: an identity operation (C1/E), inversion center (Ci), mirror plane (Cs), rotation axis (C2,C3,...), or improper rotation (S4, S6, ...) """ def __init__(self, symbol, matrix, **kwargs): self.symbol = symbol self.matrix = matrix for kw, val in kwargs.iteritems(): setattr(self, kw, val)
[docs] def get_axis(self): """ Returns normal of the plane for Cs or axis of rotation for a Cn :param symm_elem: symmetry element (with attributes 'symbol' and 'matrix') :return: array of shape (3,) """ mat = np.identity(4) mat[:3, :3] = self.matrix symbol = self.symbol if symbol == 'Cs': point, normal = trns.reflection_from_matrix(mat) assert np.allclose(point[:3], np.zeros(3)) return normal elif symbol[0] == 'C' and symbol[-1].isdigit(): angle, normal, point = trns.rotation_from_matrix(mat) assert np.allclose(point[:3], np.zeros(3)) return normal elif symbol[0] == 'S' and symbol[-1].isdigit(): normal = improper_axis_from_matrix(self.matrix) return normal else: raise ValueError('Unrecognized symmetry type %s' % self.symbol)
[docs]class MolecularSymmetry(object): def __init__(self, mol, symbol, rms, orientation=None, elems=None, **kwargs): self.mol = mol self.symbol = symbol self.rms = rms self.orientation = mdt.utils.if_not_none(orientation, mol.atoms.position) self.elems = mdt.utils.if_not_none(elems, []) for kw, val in kwargs.iteritems(): setattr(self, kw, val)
[docs] def get_symmetrized_coords(self, elem): """ Symmetrize the molecule based on the symmetry operation This will work as long as the symmetry operation brings each atom closest to a symmetry relation. """ import scipy.spatial.distance # First, apply the transformation oriented_coords = self.orientation transformed_coords = self.orientation.T.ldot(elem.matrix).T # Next, calculate the correspondence between the untransformed and transformed atoms align_to_transform = {} # map between the original positions and their transformed positions transform_to_align = {} # inverse byelement = mdt.utils.Categorizer(lambda x: x.element, self.mol.atoms) for elemname, atoms in byelement.iteritems(): indices = np.array([atom.index for atom in atoms]) atoms_aligned = oriented_coords[indices].defunits_value() atoms_transformed = transformed_coords[indices].defunits_value() distances = scipy.spatial.distance.cdist(atoms_aligned, atoms_transformed) for a_idx, t_idx in enumerate(distances.argmin(axis=0)): align_to_transform[indices[a_idx]] = indices[t_idx] for t_idx, a_idx in enumerate(distances.argmin(axis=1)): transform_to_align[indices[t_idx]] = indices[a_idx] # Make the positions exactly symmetric by averaging them pos = np.zeros(transformed_coords.shape) * u.default.length for align_atom, transform_atom in align_to_transform.iteritems(): assert transform_to_align[transform_atom] == align_atom, \ 'Molecule is too far from this symmetry to symmetrize' pos[transform_atom] = (oriented_coords[transform_atom] + transformed_coords[align_atom]) / 2.0 return pos
[docs]def improper_axis_from_matrix(matrix): """ Return rotation angle and axis / mirror plane normal from improper rotation matrix. """ R = np.array(matrix, dtype=np.float64, copy=False) R33 = R[:3, :3] # direction: unit eigenvector of R33 corresponding to eigenvalue of 1 w, W = np.linalg.eig(R33.T) i = np.where(abs(np.real(w) + 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue -1") direction = np.real(W[:, i[-1]]).squeeze() return direction