# **************************************************************************
# *
# * Authors: Roberto Marabini
# *
# * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
# *
# * This program is free software; you can redistribute it and/or modify
# * it under the terms of the GNU General Public License as published by
# * the Free Software Foundation; either version 3 of the License, or
# * (at your option) any later version.
# *
# * This program is distributed in the hope that it will be useful,
# * but WITHOUT ANY WARRANTY; without even the implied warranty of
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# * GNU General Public License for more details.
# *
# * You should have received a copy of the GNU General Public License
# * along with this program; if not, write to the Free Software
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
# * 02111-1307 USA
# *
# * All comments concerning this program package may be sent to the
# * e-mail address 'scipion@cnb.csic.es'
# *
# **************************************************************************
"""
This module returns the matrices related with the different
point symmetries. Code based on chimera file sym.py
"""
import logging
logger = logging.getLogger(__name__)
from math import sin, cos, pi, acos, sqrt, asin
import numpy as np
from numpy.linalg import inv as matrix_inverse
from numpy import dot as matrix_multiply
import operator
import pwem.constants as cts
DEBUG = False # set to True for debuging
if DEBUG:
from pwem.viewers.viewer_chimera import Chimera
# expansion factor is use to create chimera bild files
# this is the factor to multiply the coordinates
expansionFactor = 100
def _column(matrix, i):
"""Returns column i of matrix as a list."""
return [row[i] for row in matrix]
def _applyMatrix(tf, points):
"""
Args:multiply point by a matrice list
"""
tf = np.array(tf)
r = matrix_multiply(points, np.transpose(tf[:, :3]))
np.add(r, tf[:, 3], r)
return r
def __length(v):
"""Returns the length of a vector."""
d = sqrt(sum([e*e for e in v]))
return d
def _normalizeVector(v):
"""Returns a unit vector in the same direction as v."""
d = __length(v)
if d == 0:
d = 1
return tuple([e/d for e in v])
def _reflectMatrix(axis):
"""Returns a matrix that reflects points in the plane
perpendicular to axis."""
axis = _normalizeVector(axis)
return ((1 - 2*axis[0]*axis[0], -2*axis[0]*axis[1], -2*axis[0]*axis[2], 0),
(-2*axis[1]*axis[0], 1 - 2*axis[1]*axis[1], -2*axis[1]*axis[2], 0),
(-2*axis[2]*axis[0], -2*axis[2]*axis[1], 1 - 2*axis[2]*axis[2], 0))
def _rotationTransform(axis, angle, center=(0, 0, 0)):
""" returns matrix that rotates "angle" angles arround vector axis
angle is in degrees
"""
axis = _normalizeVector(axis)
arad = angle*pi/180.0
sa = sin(arad)
ca = cos(arad)
if axis == (0, 0, 1) and center == (0, 0, 0):
return ((ca, -sa, 0, 0),
(sa, ca, 0, 0),
(0, 0, 1, 0))
k = 1 - ca
ax, ay, az = axis
tf = ((1 + k*(ax*ax-1), -az*sa+k*ax*ay, ay*sa+k*ax*az, 0),
(az*sa+k*ax*ay, 1 + k*(ay*ay-1), -ax*sa+k*ay*az, 0),
(-ay*sa+k*ax*az, ax*sa+k*ay*az, 1 + k*(az*az-1), 0))
if center == (0, 0, 0):
return tf
cx, cy, cz = center
c_tf = ((1, 0, 0, cx), (0, 1, 0, cy), (0, 0, 1, cz))
inv_c_tf = ((1, 0, 0, -cx), (0, 1, 0, -cy), (0, 0, 1, -cz))
rtf = _multiplyMatrices(c_tf, tf, inv_c_tf)
return rtf
def _translationMatrix(shift):
""" returns matriz that shift point by vector shift"""
tf = np.array(((1.0, 0, 0, shift[0]),
(0, 1.0, 0, shift[1]),
(0, 0, 1.0, shift[2])))
return tf
def _identityMatrix():
""" returns identity matrix"""
return (1.0, 0, 0, 0), (0, 1.0, 0, 0), (0, 0, 1.0, 0)
def _invertMatrix(tf):
""" returns inverse 3x3 matrix"""
tf = np.array(tf)
r = tf[:, :3]
t = tf[:, 3]
tfinv = np.zeros((3, 4), np.float)
rinv = tfinv[:, :3]
tinv = tfinv[:, 3]
rinv[:, :] = matrix_inverse(r)
tinv[:] = matrix_multiply(rinv, -t)
return tfinv
def _multiplyMatrices(*mlist):
""" returns matrix product of matrices in mlist"""
if len(mlist) == 2:
m1, m2 = mlist
p = [[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0]]
for r in range(3):
for c in range(4):
p[r][c] = m1[r][0]*m2[0][c] +\
m1[r][1]*m2[1][c] +\
m1[r][2]*m2[2][c]
p[r][3] += m1[r][3]
p = tuple(map(tuple, p))
else:
p = _multiplyMatrices(*mlist[1:])
p = _multiplyMatrices(mlist[0], p)
return p
def _matrixProducts(mlist1, mlist2):
""" returns list of matrix products of matrices in mlist1 and mlist2"""
plist = []
for m1 in mlist1:
for m2 in mlist2:
m1xm2 = _multiplyMatrices(m1, m2)
plist.append(m1xm2)
return plist
def _coordinateTransformList(tflist, ctf):
""" returns list of coordinate transforms of tflist by ctf"""
ctfinv = _invertMatrix(ctf)
return [_multiplyMatrices(ctfinv, tf, ctf) for tf in tflist]
def _recenterSymmetries(tflist, center):
if center == (0, 0, 0):
return tflist
ctf = _translationMatrix([-x for x in center])
return _coordinateTransformList(tflist, ctf)
def _transposeMatrix(tf):
""" returns transpose of 3x3 matrix"""
return ((tf[0][0], tf[1][0], tf[2][0], tf[0][3]),
(tf[0][1], tf[1][1], tf[2][1], tf[1][3]),
(tf[0][2], tf[1][2], tf[2][2], tf[2][3]))
# ==================== End of utility functions ====================
[docs]class SymmetryHelper:
""" This is a static class so there is no need to instantiate it.
Given a symmetry the class provides methods to "move" projection
directions to the equivalent orientation that lays within the unit
cells defined at https://scipion-em.github.io/docs/release-3.0.0/docs/developer/symmetries/symmetries.html?highlight=symmetry."""
# Dictionary to cache all the symmetry matrices and unit cell
# planes for a particular symmetry
_matricesAndPlanes = {}
[docs] @classmethod
def getSymmetryKey(cls,sym, n):
return "%s_%s" % (sym, n)
[docs] @classmethod
def getSymmetryMatricesAndPlanes(cls, sym, n):
""" Lazy return a tuple of symmetry matrices and unit cell planes for the symmetry and order passed
:param sym: Symmetry type
:param n: symmetry order
"""
key = cls.getSymmetryKey(sym, n)
if key not in cls._matricesAndPlanes.keys():
logger.info("Gettting symmetry matrices and unitcell planes for %s, %s." % (sym,n))
matrixSet = getSymmetryMatrices(sym=sym, n=n)
# get unit cell vectors
_, unitCellPlanes = getUnitCell(sym=sym, n=n, generalize=False)
cls._matricesAndPlanes[key] = (matrixSet, unitCellPlanes)
return cls._matricesAndPlanes[key]
[docs] @classmethod
def moveParticleInsideUnitCell(cls, particle, symmetry, symmetryOrder):
""" apply to the projection direction the symmetry matrix that "moves" the projection
direction to the equivalent orientation that lays within the "canonical" unit
This is useful to move for example Cryosparc angles to a same unit cell
and therefore been able to plot angular distribution. This method can be used in the typical "updateParticle"
callback when generating any output set.
:param particle: Particle to move.
:param symmetry: symmetry type in scipion's convention.
:param symmetryOrder: order of the symmetry
"""
# For C1 particles we do not move angles to unit cell.
if symmetry==cts.SYM_CYCLIC and symmetryOrder==1:
logger.debug("Cancelling unit cell migration due to C1 symmetry.")
return
matrices, unitCellPlanes = cls.getSymmetryMatricesAndPlanes(symmetry, symmetryOrder)
cls._moveParticleInsideUnitCell(particle, matrices, unitCellPlanes)
@classmethod
def _moveParticleInsideUnitCell(cls, particle, matrixSet, unitCellPlanes):
""" Move the particle projection direction so it lays inside the canonical unit cell,
symmetry matrices and unit cell planes needs to be precomputed
:param particle: to move
:param matrixSet: value returned when calling getSymmetryMatrices with the right symmetry and order
:param unitCellPlanes: second term returned by getUnitCell when called with the right symmetry and order
"""
# get projection direction
colunm = particle.getTransform().getMatrix()[0:3, 2]
u, v, w = unitCellPlanes
if (np.dot(colunm, u) > 0) and \
(np.dot(colunm, v) > 0) and \
(np.dot(colunm, w) > 0):
# particle is inside unit cell so nothing needs to be done
return particle
else:
# loop over symmetry matrices
for i, matrix in enumerate(matrixSet):
matrix3 = np.array(matrix)
matrix3 = np.delete(matrix3, -1, 1)
columPrime = matrix3.dot(colunm)[:3]
if (np.dot(columPrime, u) > 0) and \
(np.dot(columPrime, v) > 0) and \
(np.dot(columPrime, w) > 0):
matrix = np.dot(matrix, particle.getTransform().getMatrix())
t = particle.getTransform()
t.setMatrix(matrix)
particle.setTransform(t)
return particle
logger.info("Error: something went wrong in moveParticlesInsideUnitCell."
" No matrix found to move the particle projection direction inside the unit cell."
" particle id: %s" % particle.getObjId())
[docs]def moveParticlesInsideUnitCell(setIN, setOUT, sym=cts.SYM_CYCLIC, n=1):
""" Move particles projection direction inside the unit cell associated
to symmetry sym.
This function (1) gets the symmetry matrices and (2) applies them to
each projection direction until is inside the unit cell.
"""
counter = 0 # counter for the number of particles
totalNumberOfParticles = len(setIN)
# for each particle in set
for par in setIN:
counter += 1
if counter % 1000 == 0:
# print something to keep the user engaged
logger.info(f"{counter}/{totalNumberOfParticles}")
particle = SymmetryHelper.moveParticleInsideUnitCell(par, sym, n)
if particle:
setOUT.append(particle)
[docs]def moveParticleInsideUnitCell(particle, matrixSet, unitCellPlanes):
logger.warning("FOR DEVELOPERS: This method is deprecated. Call SymmetryHelper.moveParticleInsideUnitCell for convenience.")
SymmetryHelper._moveParticleInsideUnitCell(particle, matrixSet, unitCellPlanes)
[docs]def getSymmetryMatrices(sym=cts.SYM_CYCLIC,
n=1,
circumscribed_radius=1,
center=(0, 0, 0),
offset=None):
""" interface between scipion and chimera code
chimera code uses tuples of tuples as matrices
but scipion uses np.arrays (lists of lists)
so let us convert them here.
Direct use of the classes return 3x3 arrays
which may be OK in many cases, but not in all
"""
if sym == cts.SYM_CYCLIC:
c = Cyclic(n=n,
center=center,
circumscribed_radius=circumscribed_radius,
offset=offset)
matrices = c.symmetryMatrices()
elif sym == cts.SYM_DIHEDRAL_X or sym == cts.SYM_DIHEDRAL_Y:
d = Dihedral(n=n,
center=center,
circumscribed_radius=circumscribed_radius,
offset=offset, sym=sym)
matrices = d.symmetryMatrices()
elif sym == cts.SYM_OCTAHEDRAL:
o = Octahedral(center=center,
circumscribed_radius=circumscribed_radius,
offset=offset)
matrices = o.symmetryMatrices()
elif (sym == cts.SYM_TETRAHEDRAL_222 or
sym == cts.SYM_TETRAHEDRAL_Z3 or
sym == cts.SYM_TETRAHEDRAL_Z3R):
t = Tetrahedral(center=center,
circumscribed_radius=circumscribed_radius,
sym=sym)
matrices = t.symmetryMatrices()
elif (sym == cts.SYM_I222 or sym == cts.SYM_I222r or
sym == cts.SYM_In25 or sym == cts.SYM_In25r or
sym == cts.SYM_I2n3 or sym == cts.SYM_I2n3r or
sym == cts.SYM_I2n5 or sym == cts.SYM_I2n5r):
matrices = __icosahedralSymmetryMatrices(sym, center)
# '222' 2-fold symmetry along x, y, and z axes.
# '222r' '222' with 90 degree rotation around z.
# '2n5' 2-fold symmetry along x and 5-fold along z.
# '2n5r' '2n5' with 180 degree rotation about y.
# 'n25' 2-fold symmetry along y and 5-fold along z.
# 'n25r' 'n25' with 180 degree rotation about x.
# '2n3' 2-fold symmetry along x and 3-fold along z.
# '2n3r' '2n3' with 180 degree rotation about y.
# convert from 4x3 to 4x4 matrix, Scipion standard
extraRow = (0., 0., 0., 1.)
for i in range(len(matrices)):
matrices[i] += (extraRow,)
# convert from sets to lists Scipion standard
return np.array(matrices)
[docs]def getUnitCell(sym=cts.SYM_CYCLIC,
n=1,
circumscribed_radius=1,
center=(0, 0, 0),
offset=None,
generalize=False):
""" return vectors normal to the unit cell faces
"""
if sym == cts.SYM_CYCLIC:
c = Cyclic(n=n,
center=center,
circumscribed_radius=circumscribed_radius,
offset=offset)
vectorsEdge, vectorsPlane = c.unitCellPlanes()
# This function return a three vector planes, but cyclic symmetry only
# need two, so we duplicate the last one
vectorsPlane.append(vectorsPlane[-1])
if DEBUG:
bildFileName = f'/tmp/C{n}.bild'
vLabels = ['v1', 'v2', 'v3']
pLabels = ['p1', 'p2']
elif sym == cts.SYM_DIHEDRAL_X or sym == cts.SYM_DIHEDRAL_Y:
logger.info("getUnitCell: %s" % cts.SCIPION_SYM_NAME[sym])
d = Dihedral(n=n,
center=center,
circumscribed_radius=circumscribed_radius,
offset=offset, sym=sym)
vectorsEdge, vectorsPlane = d.unitCellPlanes()
if DEBUG:
vLabels = ['v1', 'v2', 'eigenvector']
pLabels = ['p1', 'p2', 'p3']
elif sym == cts.SYM_OCTAHEDRAL:
o = Octahedral(center=center,
circumscribed_radius=circumscribed_radius)
vectorsEdge, vectorsPlane = o.unitCellPlanes()
if DEBUG:
vLabels = ['v1', 'v2', 'v3']
pLabels = ['p1', 'p2', 'p3']
elif (sym == cts.SYM_TETRAHEDRAL_Z3 or
sym == cts.SYM_TETRAHEDRAL or
sym == cts.SYM_TETRAHEDRAL_Z3R):
t = Tetrahedral(center=center,
circumscribed_radius=circumscribed_radius,
sym=sym)
vectorsEdge, vectorsPlane = t.unitCellPlanes()
if DEBUG:
vLabels = ['v1', 'v2', 'v3']
pLabels = ['p1', 'p2', 'p3']
elif (sym == cts.SYM_I222 or sym == cts.SYM_I222r or
sym == cts.SYM_In25 or sym == cts.SYM_In25r or
sym == cts.SYM_I2n3 or sym == cts.SYM_I2n3r or
sym == cts.SYM_I2n5 or sym == cts.SYM_I2n5r):
i_sym = cts.SCIPION_SYM_NAME[sym][1:]
i = Icosahedron(center=center,
circumscribed_radius=circumscribed_radius,
orientation=i_sym)
vectorsEdge, vectorsPlane = i.unitCellPlanes()
if DEBUG:
vLabels = ['v1', 'v2', 'v3']
pLabels = ['p1', 'p2', 'p3']
if DEBUG:
bildFileName = cts.SCIPION_SYM_NAME[sym] + '.bild'
Chimera.createCoordinateAxisFile(dim=expansionFactor,
bildFileName=bildFileName,
sampling=1)
colors = ['cyan', 'green', 'magenta', 'navy blue']
with open(bildFileName, 'a') as f:
for i, vector in enumerate(vectorsEdge):
f.write(f'.comment {vLabels[i]}\n')
f.write(f'.color {colors[i]}\n')
x = vector[0]*expansionFactor
y = vector[1]*expansionFactor
z = vector[2]*expansionFactor
f.write(f'.arrow 0 0 0 {x} {y} {z} 0.20 0.40 0.75\n')
for i, vector in enumerate(vectorsPlane):
f.write(f'.comment {pLabels[i]}\n')
f.write(f'.color {colors[i]}\n')
x = vector[0]
y = vector[1]
z = vector[2]
r = expansionFactor/2
f.write(f'.cone 0 0 0 {x/(expansionFactor*100)} \
{y/(expansionFactor*100)} \
{z/(expansionFactor*100)} {r}\n')
f.write(f'.arrow 0 0 0 {x*expansionFactor} \
{y*expansionFactor} \
{z*expansionFactor} 0.20 0.40 0.75\n')
f.close()
if generalize:
for i, v in enumerate(vectorsEdge):
vectorsEdge[i] = np.append(v, 1.)
for i, v in enumerate(vectorsPlane):
vectorsPlane[i] = np.append(v, 1.)
return vectorsEdge, vectorsPlane
def __icosahedralSymmetryMatrices(orientation=cts.SYM_I222, center=(0, 0, 0)):
if orientation == cts.SYM_I222:
sym = '222'
elif orientation == cts.SYM_I222r:
sym = '222r'
elif orientation == cts.SYM_In25:
sym = 'n25'
elif orientation == cts.SYM_In25r:
sym = 'n25r'
elif orientation == cts.SYM_I2n5:
sym = '2n5'
elif orientation == cts.SYM_I2n5r:
sym = '2n5r'
elif orientation == cts.SYM_I2n3:
sym = '2n3'
elif orientation == cts.SYM_I2n3r:
sym = '2n3r'
i = Icosahedron(orientation=sym, center=center)
return list(i.icosahedralSymmetryMatrices())
def __icosahedralUnitCellPlanes(orientation=cts.SYM_I222, center=(0, 0, 0)):
if orientation == cts.SYM_I222:
sym = '222'
elif orientation == cts.SYM_I222r:
sym = '222r'
elif orientation == cts.SYM_In25:
sym = 'n25'
elif orientation == cts.SYM_In25r:
sym = 'n25r'
elif orientation == cts.SYM_I2n5:
sym = '2n5'
elif orientation == cts.SYM_I2n5r:
sym = '2n5r'
elif orientation == cts.SYM_I2n3:
sym = '2n3'
elif orientation == cts.SYM_I2n3r:
sym = '2n3r'
i = Icosahedron(orientation=sym, center=center)
return list(i.unitCellPlanes())
[docs]class Cyclic(object):
"""cyclic class.
Allow to compute symmetry matrices, vertices,
symmetry axis and unit cell planes
"""
def __init__(self, n, circumscribed_radius=1, center=(
0, 0, 0), offset=None):
"""
:Parameters:
n: int order of symmetry
circumscribed_radius: float radius of the circumscribed sphere
center: tuple of 3 floats, center of the circumscribed sphere
offset: float, angle in degrees to rotate the symmetry axis
it modifies the unit cell but not the symmetry matrices
"""
self.n = n
self.circumscribed_radius = circumscribed_radius
self.center = center
self.matrices = None
self.offset = offset
[docs] def symmetryMatrices(self):
""" get Matrices for cyclic symmetry of order n
"""
if self.matrices is not None:
return self.matrices
self.matrices = []
for k in range(self.n):
a = 2*pi * np.float32(k) / self.n
c = cos(a)
s = sin(a)
tf = ((c, -s, 0, 0),
(s, c, 0, 0),
(0, 0, 1, 0))
self.matrices.append(tf)
self.matrices = _recenterSymmetries(self.matrices, self.center)
return self.matrices
[docs] def unitCellPlanes(self):
""" get planes that define a unit cell for cyclic symmetry of order n
"""
if self.offset is None:
self.offset = pi/(self.n)
_ = self.symmetryMatrices()
# these three vectors are enges of the unit cell
angle = -2*pi/self.n
a1 = self.offset
a2 = angle + self.offset
c1 = cos(a1)
s1 = sin(a1)
c2 = cos(a2)
s2 = sin(a2)
v1 = np.array([c1, s1, 0.])
v2 = np.array([c2, s2, 0.])
v3 = np.array([0., 0., 1.]) # keep this dim 3 instead of dim 4
# cross product of v1/v2 with eigenvetor
# these two vectors are normal to the planes that define the unit cell
plane1 = np.cross(v1, v3)
plane2 = np.cross(v3, v2)
if DEBUG:
logger.info("v1 %s" % v1)
logger.info("v2 %s" % v2)
logger.info("v3 %s" % v3)
logger.info("plane1 %s" % plane1)
logger.info("plane2 %s" % plane2)
c = self.circumscribed_radius
return [v1*c, v2*c, v3*c], [plane1*c, plane2*c]
[docs]class Dihedral(Cyclic):
"""dihedral class.
Allow to compute symmetry matrices, vertices,
symmetry axis and unit cell planes
"""
def __init__(self, n, circumscribed_radius=1, center=(
0, 0, 0), offset=None, sym=cts.SYM_DIHEDRAL_X):
"""
:Parameters:
n: int order of symmetry
circumscribed_radius: float radius of the circumscribed sphere
center: tuple of 3 floats, center of the circumscribed sphere
offset: float, angle in degrees to rotate the symmetry axis
it modifies the unit cell but not the symmetry matrices
"""
super().__init__(n=n,
circumscribed_radius=circumscribed_radius,
center=center,
offset=offset)
self.sym = sym
[docs] def symmetryMatrices(self):
clist = super().symmetryMatrices()
if self.sym == cts.SYM_DIHEDRAL_X:
# matrix that reflect in the x axis
reflect = ((1, 0, 0, 0), (0, -1, 0, 0), (0, 0, -1, 0))
else:
logger.info("Be carefull untested code")
# matrix that reflect in the y axis
reflect = ((-1, 0, 0, 0), (0, 1, 0, 0), (0, 0, -1, 0))
tflist = _matrixProducts([_identityMatrix(), reflect], clist)
tflist = _recenterSymmetries(tflist, self.center)
return tflist
[docs] def unitCellPlanes(self):
""" get planes that define a unit cell for dihedral symmetry of order n
"""
if self.offset is None:
self.offset = pi/(self.n) + pi/2
vectorsEdge, vectorsPlane = super().unitCellPlanes()
vectorsPlane.append(np.array([0., 0., self.circumscribed_radius]))
if DEBUG:
logger.info("v1 %s" % vectorsEdge[0])
logger.info("v2 %s" % vectorsEdge[1])
logger.info("eigenvector %s" % vectorsEdge[2])
logger.info("plane1 %s" % vectorsPlane[0])
logger.info("plane2 %s" % vectorsPlane[1])
logger.info("plane3 %s" % vectorsPlane[2])
return vectorsEdge, vectorsPlane
[docs]class Tetrahedral(object):
"""Tetrahedral class.
Allow to compute symmetry matrices, vertices,
symmetry axis and unit cell planes
"""
def __init__(self, circumscribed_radius=1, center=(
0, 0, 0), sym=cts.SYM_TETRAHEDRAL_Z3):
"""
:Parameters:
n: int order of symmetry
circumscribed_radius: float radius of the circumscribed sphere
center: tuple of 3 floats, center of the circumscribed sphere
"""
self.circumscribed_radius = circumscribed_radius
self.center = center
self.matrices = None
self.sym = sym
[docs] def symmetryMatrices(self):
"""
identity
4 * rotation by 120 clockwise (seen from a vertex):
(234), (143), (412), (321)
4 * rotation by 120 counterclockwise (ditto)
3 * rotation by 180
"""
# simmetries are rotations around theses axis by these angles
self.aa = (((0, 0, 1), 0), ((1, 0, 0), 180), # 0 1
((0, 1, 0), 180), ((0, 0, 1), 180), # 2 3
((1, 1, 1), 120), ((1, 1, 1), 240), # 4 5
((-1, -1, 1), 120), ((-1, -1, 1), 240), # 6 7
((-1, 1, -1), 120), ((-1, 1, -1), 240), # 8 9
((1, -1, -1), 120), ((1, -1, -1), 240)) # 10 11
syms = [_rotationTransform(axis, angle) for axis, angle in self.aa]
if self.sym == cts.SYM_TETRAHEDRAL_Z3R:
# convention, 3-fold on z, 3-fold in yz plane sign(y) = sign(z)
self.tf = _multiplyMatrices(
_rotationTransform((0, 0, 1), -45.0),
_rotationTransform((1, 0, 0),
-acos(1 / sqrt(3)) * 180 / pi))
syms = _coordinateTransformList(syms, self.tf)
elif self.sym == cts.SYM_TETRAHEDRAL_Z3:
# convention, 3-fold on z, 3-fold in yz plane sign(y) != sign(z)
tfaux = _multiplyMatrices(
_rotationTransform((0, 0, 1), -45.0),
_rotationTransform((1, 0, 0),
-acos(1 / sqrt(3)) * 180 / pi)
)
self.tf = _multiplyMatrices(
tfaux,
_reflectMatrix((0, 0, 1))
)
syms = _coordinateTransformList(syms, self.tf)
syms = _recenterSymmetries(syms, self.center)
return syms
[docs] def unitCellPlanes(self):
""" get planes that define a unit cell for tetrahedral symmetry of
order n
"""
logger.info("Unit cell sym %s" % cts.SCIPION_SYM_NAME[self.sym])
_ = self.symmetryMatrices()
# four corners tetahedron
v0 = _normalizeVector(np.array(self.aa[4][0])) *\
self.circumscribed_radius
v1 = _normalizeVector(np.array(self.aa[6][0])) *\
self.circumscribed_radius
v2 = _normalizeVector(np.array(self.aa[8][0])) *\
self.circumscribed_radius
v3 = _normalizeVector(np.array(self.aa[10][0])) *\
self.circumscribed_radius
if self.sym == cts.SYM_TETRAHEDRAL_Z3 or\
self.sym == cts.SYM_TETRAHEDRAL_Z3R:
self.tf = _invertMatrix(self.tf)
if self.sym == cts.SYM_TETRAHEDRAL_Z3:
_3fold_1 = np.dot(self.tf, np.append(v0, 1))
_3fold_2 = np.dot(self.tf, np.append(v1, 1))
_3fold_3 = np.dot(self.tf, np.append(v3, 1))
else:
_3fold_1 = np.dot(self.tf, np.append(v1, 1))
_3fold_2 = np.dot(self.tf, np.append(v0, 1))
_3fold_3 = np.dot(self.tf, np.append(v3, 1))
else:
_3fold_1 = v1
_3fold_2 = v0
_3fold_3 = v2
# colors = ['cyan', 'green', 'magenta', 'navy blue']
vp = np.add(_3fold_1, _3fold_2)
vp = _normalizeVector(np.add(vp, _3fold_3))
plane1 = _normalizeVector(np.cross(_3fold_1, _3fold_2)) # cyan
plane2 = _normalizeVector(np.cross(_3fold_2, vp)) # green
plane3 = _normalizeVector(np.cross(vp, _3fold_1)) # magenta
if DEBUG:
logger.info("v0 %s" % _3fold_1)
logger.info("v1 %s" % _3fold_2)
logger.info("v2 %s" % vp)
logger.info("plane1 %s" % plane1)
logger.info("plane2 %s" % plane2)
logger.info("plane3 %s" % plane3)
return [_3fold_1, vp, _3fold_2], [plane1, plane2, plane3]
[docs]class Octahedral(object):
"""Octahedral class.
Allow to compute symmetry matrices,
symmetry axis and unit cell planes of an octahedral symmetry.
4-folds along x, y, z axes.
"""
def __init__(self, circumscribed_radius=1, center=(
0, 0, 0), offset=None):
"""
:Parameters:
circumscribed_radius: float radius of the circumscribed sphere
center: tuple of 3 floats, center of the circumscribed sphere
offset: float, angle in degrees to rotate the symmetry axis
it modifies the unit cell but not the symmetry matrices
"""
self.circumscribed_radius = circumscribed_radius
self.center = center
self.matrices = None
self.offset = offset
[docs] def symmetryMatrices(self):
""" get Matrices for cyclic symmetry of order n
"""
self.c4 = (((0, 0, 1), 0),
((0, 0, 1), 90),
((0, 0, 1), 180),
((0, 0, 1), 270))
self.cube = (((1, 0, 0), 0),
((1, 0, 0), 90),
((1, 0, 0), 180),
((1, 0, 0), 270),
((0, 1, 0), 90),
((0, 1, 0), 270))
c4syms = [_rotationTransform(axis, angle)
for axis, angle in self.c4]
cubesyms = [_rotationTransform(axis, angle)
for axis, angle in self.cube]
syms = _matrixProducts(cubesyms, c4syms)
self.matrices = _recenterSymmetries(syms, self.center)
return self.matrices
[docs] def unitCellPlanes(self):
""" get planes that define a unit cell for an octahedral symmetry
"""
_ = self.symmetryMatrices()
# four corners tetahedron
_4fold_1 = _normalizeVector(
np.array([0., 0., 1.])) * self.circumscribed_radius
_3fold_1 = _normalizeVector(
np.array([1., 1., 1.])) * self.circumscribed_radius
_3fold_2 = _normalizeVector(
np.array([-1., 1., 1.])) * self.circumscribed_radius
plane1 = _normalizeVector(
np.cross(_3fold_2, _4fold_1)) # cyan
plane2 = _normalizeVector(
np.cross(_3fold_1, _3fold_2)) # green
plane3 = _normalizeVector(
np.cross(_4fold_1, _3fold_1)) # magenta
if DEBUG:
logger.info("v0 %s" % _4fold_1)
logger.info("v1 %s" % _3fold_2)
logger.info("v2 %s" % _3fold_1)
logger.info("plane1 %s" % plane1)
logger.info("plane2 %s" % plane2)
logger.info("plane3 %s" % plane3)
return [_4fold_1, _3fold_2, _3fold_1], [plane1, plane2, plane3]
# TODO: Why is this a global variable instead of
# a class variable?
icos_matrices = {} # Maps orientation name to 60 matrices.
[docs]class Icosahedron(object):
"""Icosahedron class.
Allow to compute symmetry matrices, icosahedron vertices,
symmetry axis and unit cell planes
"""
def __init__(self, circumscribed_radius=1, orientation='222', center=(
0, 0, 0)):
"""point = np.array([0, 1, PHI]"""
self.circumscribed_radius = circumscribed_radius
self.orientation = orientation
self.center = center
# Triangle edge length of unit icosahedron.
self.e = sqrt(2 - 2 / sqrt(5))
self.vertices = None # icosahedron vertices
self.triangles = None # icosahedron faces
self._3foldAxis = []
self._2foldAxis = []
self.edges = None # icosahedron edges
# ---------------------------------------------------------------------------------------------
# Coordinates systems.
# '222' 2-fold symmetry along x, y, and z axes.
# '222r' '222' with 90 degree rotation around z.
# '2n5' 2-fold symmetry along x and 5-fold along z.
# '2n5r' '2n5' with 180 degree rotation about y.
# 'n25' 2-fold symmetry along y and 5-fold along z.
# 'n25r' 'n25' with 180 degree rotation about x.
# '2n3' 2-fold symmetry along x and 3-fold along z.
# '2n3r' '2n3' with 180 degree rotation about y.
#
self.coordinate_system_names = ('222', '222r', '2n5',
'2n5r', 'n25', 'n25r', '2n3', '2n3r')
self.icosahedronEdgeLength = \
self.icosahedronEdgeLength(circumscribed_radius)
self.angle23, self.angle25, self.angle35 = self.icosahedronAngles()
self.icosahedronGeometry()
# ---------------------------------------------------------------------------------------------
# 60 icosahedral symmetry matrices.
#
[docs] def icosahedralSymmetryMatrices(self):
t = self.icosahedralMatrixTable()
tflist = _recenterSymmetries(t.get(self.orientation, None),
self.center)
return tflist
# -------------------------------------------------------------------------
# Edge length of icosahedron with a certain radio of the circumscribed
# sphere:
# According to Radio of the circumscribed sphere = 1 (vertices),
# the icosahedronEdgeLength should be 1.0515
[docs] def icosahedronEdgeLength(self, circumscribed_radius):
return 4 * circumscribed_radius / sqrt(10 + 2 * sqrt(5))
# -------------------------------------------------------------------------
# Matrices for mapping between different icosahedron coordinate frames.
#
# ---------------------------------------------------------------------------------------
# Compute icosahedral transformation matrices for different
# coordinate systems.
#
[docs] def icosahedralMatrixTable(self):
global icos_matrices
if icos_matrices:
return icos_matrices
c = cos(2 * pi / 5) # .309016994
c2 = cos(4 * pi / 5) # -.809016994
icos_matrices['222'] = (
((1.0, 0.0, 0.0, 0.0),
(0.0, 1.0, 0.0, 0.0),
(0.0, 0.0, 1.0, 0.0)),
((c2, -0.5, c, 0.0),
(-0.5, c, c2, 0.0),
(c, c2, -0.5, 0.0)),
((0.0, 1.0, 0.0, 0.0),
(0.0, 0.0, -1.0, 0.0),
(-1.0, 0.0, 0.0, 0.0)),
((-c2, -0.5, -c, 0.0),
(-0.5, -c, c2, 0.0),
(c, -c2, -0.5, 0.0)),
((0.5, c, c2, 0.0),
(-c, c2, -0.5, 0.0),
(c2, 0.5, -c, 0.0)),
((-c, c2, -0.5, 0.0),
(c2, 0.5, -c, 0.0),
(0.5, c, c2, 0.0)),
((c2, 0.5, -c, 0.0),
(0.5, c, c2, 0.0),
(-c, c2, -0.5, 0.0)),
((c2, -0.5, -c, 0.0),
(0.5, -c, c2, 0.0),
(c, c2, 0.5, 0.0)),
((-c, -c2, -0.5, 0.0),
(c2, -0.5, -c, 0.0),
(-0.5, c, -c2, 0.0)),
((0.5, -c, c2, 0.0),
(-c, -c2, -0.5, 0.0),
(-c2, 0.5, c, 0.0)),
((0.0, 0.0, -1.0, 0.0),
(-1.0, 0.0, 0.0, 0.0),
(0.0, 1.0, 0.0, 0.0)),
((-0.5, -c, c2, 0.0),
(c, -c2, -0.5, 0.0),
(-c2, -0.5, -c, 0.0)),
((-0.5, c, c2, 0.0),
(c, c2, -0.5, 0.0),
(c2, -0.5, c, 0.0)),
((-c, c2, -0.5, 0.0),
(-c2, -0.5, c, 0.0),
(-0.5, -c, -c2, 0.0)),
((c2, 0.5, -c, 0.0),
(-0.5, -c, -c2, 0.0),
(c, -c2, 0.5, 0.0)),
((0.5, c, c2, 0.0),
(c, -c2, 0.5, 0.0),
(-c2, -0.5, c, 0.0)),
((-0.5, c, c2, 0.0),
(-c, -c2, 0.5, 0.0),
(-c2, 0.5, -c, 0.0)),
((0.0, 0.0, -1.0, 0.0),
(1.0, 0.0, 0.0, 0.0),
(0.0, -1.0, 0.0, 0.0)),
((-0.5, -c, c2, 0.0),
(-c, c2, 0.5, 0.0),
(c2, 0.5, c, 0.0)),
((0.0, -1.0, 0.0, 0.0),
(0.0, 0.0, 1.0, 0.0),
(-1.0, 0.0, 0.0, 0.0)),
((c2, 0.5, c, 0.0),
(0.5, c, -c2, 0.0),
(c, -c2, -0.5, 0.0)),
((-c2, 0.5, -c, 0.0),
(0.5, -c, -c2, 0.0),
(c, c2, -0.5, 0.0)),
((-c, -c2, -0.5, 0.0),
(-c2, 0.5, c, 0.0),
(0.5, -c, c2, 0.0)),
((0.5, -c, c2, 0.0),
(c, c2, 0.5, 0.0),
(c2, -0.5, -c, 0.0)),
((c2, -0.5, -c, 0.0),
(-0.5, c, -c2, 0.0),
(-c, -c2, -0.5, 0.0)),
((-c, c2, 0.5, 0.0),
(c2, 0.5, c, 0.0),
(-0.5, -c, c2, 0.0)),
((-c, -c2, 0.5, 0.0),
(-c2, 0.5, -c, 0.0),
(-0.5, c, c2, 0.0)),
((1.0, 0.0, 0.0, 0.0),
(0.0, -1.0, 0.0, 0.0),
(0.0, 0.0, -1.0, 0.0)),
((c, -c2, -0.5, 0.0),
(-c2, -0.5, -c, 0.0),
(-0.5, -c, c2, 0.0)),
((c, c2, -0.5, 0.0),
(c2, -0.5, c, 0.0),
(-0.5, c, c2, 0.0)),
((-1.0, 0.0, 0.0, 0.0),
(0.0, 1.0, 0.0, 0.0),
(0.0, 0.0, -1.0, 0.0)),
((-c2, 0.5, -c, 0.0),
(-0.5, c, c2, 0.0),
(-c, -c2, 0.5, 0.0)),
((0.0, -1.0, 0.0, 0.0),
(0.0, 0.0, -1.0, 0.0),
(1.0, 0.0, 0.0, 0.0)),
((c2, 0.5, c, 0.0),
(-0.5, -c, c2, 0.0),
(-c, c2, 0.5, 0.0)),
((-0.5, -c, -c2, 0.0),
(-c, c2, -0.5, 0.0),
(-c2, -0.5, c, 0.0)),
((c, -c2, 0.5, 0.0),
(c2, 0.5, -c, 0.0),
(-0.5, -c, -c2, 0.0)),
((-c2, -0.5, c, 0.0),
(0.5, c, c2, 0.0),
(c, -c2, 0.5, 0.0)),
((-c2, 0.5, c, 0.0),
(0.5, -c, c2, 0.0),
(-c, -c2, -0.5, 0.0)),
((c, c2, 0.5, 0.0),
(c2, -0.5, -c, 0.0),
(0.5, -c, c2, 0.0)),
((-0.5, c, -c2, 0.0),
(-c, -c2, -0.5, 0.0),
(c2, -0.5, -c, 0.0)),
((0.0, 0.0, 1.0, 0.0),
(-1.0, 0.0, 0.0, 0.0),
(0.0, -1.0, 0.0, 0.0)),
((0.5, c, -c2, 0.0),
(c, -c2, -0.5, 0.0),
(c2, 0.5, c, 0.0)),
((0.5, -c, -c2, 0.0),
(c, c2, -0.5, 0.0),
(-c2, 0.5, -c, 0.0)),
((c, -c2, 0.5, 0.0),
(-c2, -0.5, c, 0.0),
(0.5, c, c2, 0.0)),
((-c2, -0.5, c, 0.0),
(-0.5, -c, -c2, 0.0),
(-c, c2, -0.5, 0.0)),
((-0.5, -c, -c2, 0.0),
(c, -c2, 0.5, 0.0),
(c2, 0.5, -c, 0.0)),
((0.5, -c, -c2, 0.0),
(-c, -c2, 0.5, 0.0),
(c2, -0.5, c, 0.0)),
((0.0, 0.0, 1.0, 0.0),
(1.0, 0.0, 0.0, 0.0),
(0.0, 1.0, 0.0, 0.0)),
((0.5, c, -c2, 0.0),
(-c, c2, 0.5, 0.0),
(-c2, -0.5, -c, 0.0)),
((0.0, 1.0, 0.0, 0.0),
(0.0, 0.0, 1.0, 0.0),
(1.0, 0.0, 0.0, 0.0)),
((-c2, -0.5, -c, 0.0),
(0.5, c, -c2, 0.0),
(-c, c2, 0.5, 0.0)),
((c2, -0.5, c, 0.0),
(0.5, -c, -c2, 0.0),
(-c, -c2, 0.5, 0.0)),
((c, c2, 0.5, 0.0),
(-c2, 0.5, c, 0.0),
(-0.5, c, -c2, 0.0)),
((-0.5, c, -c2, 0.0),
(c, c2, 0.5, 0.0),
(-c2, 0.5, c, 0.0)),
((-c2, 0.5, c, 0.0),
(-0.5, c, -c2, 0.0),
(c, c2, 0.5, 0.0)),
((c, -c2, -0.5, 0.0),
(c2, 0.5, c, 0.0),
(0.5, c, -c2, 0.0)),
((c, c2, -0.5, 0.0),
(-c2, 0.5, -c, 0.0),
(0.5, -c, -c2, 0.0)),
((-1.0, 0.0, 0.0, 0.0),
(0.0, -1.0, 0.0, 0.0),
(0.0, 0.0, 1.0, 0.0)),
((-c, c2, 0.5, 0.0),
(-c2, -0.5, -c, 0.0),
(0.5, c, -c2, 0.0)),
((-c, -c2, 0.5, 0.0),
(c2, -0.5, c, 0.0),
(0.5, -c, -c2, 0.0)),
)
for cs in self.coordinate_system_names:
if cs != '222':
t = self.coordinateSystemTransform(cs, '222')
tinv = self.coordinateSystemTransform('222', cs)
icos_matrices[cs] = [_multiplyMatrices(tinv, m, t)
for m in icos_matrices['222']]
return icos_matrices
# -----------------------------------------------------------------------------
[docs] def icosahedronAngles(self):
# Angles between 2-fold, 3-fold and 5-fold
# symmetry axes of an icosahedron.
angle25 = asin(self.e/2)
angle35 = asin(self.e/sqrt(3))
angle23 = asin(self.e/(2*sqrt(3)))
return angle23, angle25, angle35
[docs] def icosahedronGeometry(self):
a = 2 * self.angle25 # Angle spanned by edge from center
# 5-fold symmetry axis along z
c5 = cos(2 * pi / 5)
s5 = sin(2 * pi / 5)
tf5 = ((c5, -s5, 0, 0),
(s5, c5, 0, 0),
(0, 0, 1, 0))
# 2-fold symmetry axis along x
tf2 = ((1, 0, 0, 0),
(0, -1, 0, 0),
(0, 0, -1, 0))
p = tuple(map(operator.add, self.center,
(0, 0, self.circumscribed_radius)))
p50 = tuple(map(operator.add, self.center,
(0, self.circumscribed_radius * sin(a),
self.circumscribed_radius * cos(a))))
p51 = _applyMatrix(tf5, p50)
p52 = _applyMatrix(tf5, p51)
p53 = _applyMatrix(tf5, p52)
p54 = _applyMatrix(tf5, p53)
self.vertices = [p, p50, p51, p52, p53, p54]
self.vertices.extend([_applyMatrix(tf2, q) for q in self.vertices])
if self.orientation != '2n5':
self.tf = self.coordinateSystemTransform('2n5', self.orientation)
self.vertices = [_applyMatrix(self.tf, p) for p in self.vertices]
#
# Vertex numbering
#
# Top 1 Bottom
# 2 5 9 10
# 0 6
# 3 4 8 11
# 7
# 20 triangles composing icosahedron.
#
self.triangles = ((0, 1, 2), (0, 2, 3),
(0, 3, 4), (0, 4, 5),
(0, 5, 1), (6, 7, 8),
(6, 8, 9), (6, 9, 10),
(6, 10, 11), (6, 11, 7),
(1, 9, 2), (2, 9, 8),
(2, 8, 3), (3, 8, 7),
(3, 7, 4), (4, 7, 11),
(4, 11, 5), (5, 11, 10),
(5, 10, 1), (1, 10, 9))
for triangle in self.triangles:
self._3foldAxis.append([(item1 + item2 + item3) / 3.
for item1, item2, item3 in zip(
self.vertices[triangle[0]],
self.vertices[triangle[1]],
self.vertices[triangle[2]])])
self.edges = ((0, 1), (0, 2), (0, 3), (0, 4), (0, 5),
(1, 2), (2, 3), (3, 4), (4, 5), (5, 1),
(6, 7), (6, 8), (6, 9), (6, 10), (6, 11),
(1, 9), (2, 8), (2, 9), (3, 7), (3, 8),
(7, 8), (8, 9), (9, 10), (10, 11), (11, 7),
(4, 7), (4, 11), (5, 10), (5, 11), (1, 10))
for edge in self.edges:
self._2foldAxis.append([(item1 + item2) / 2.
for item1, item2 in zip(
self.vertices[edge[0]],
self.vertices[edge[1]])
])
[docs] def getVertices(self):
return self.vertices
[docs] def getTriangles(self):
return self.triangles
[docs] def getEdges(self):
return self.edges
[docs] def get3foldAxis(self):
return self._3foldAxis
[docs] def get2foldAxis(self):
return self._2foldAxis
[docs] def unitCellPlanes(self):
""" get planes that define a unit cell for an octahedral symmetry
"""
if self.orientation == '222':
# Unit cell in triangle 11 -> (2,8,9)
_3f = 11; _5f1 = 2; _5f2 = 8 # noqa
elif self.orientation == '222r':
# Unit cell in triangle 0 -> (0,1,2)
_3f = 0; _5f1 = 1; _5f2 = 0 # noqa
elif self.orientation == 'n25':
# Unit cell in triangle 9 -> (6,11,7)
_3f = 9; _5f1 = 6; _5f2 = 7 # noqa
elif self.orientation == 'n25r':
# Unit cell in triangle 14 -> (3,7,4)
_3f = 14; _5f1 = 3; _5f2 = 4 # noqa
elif self.orientation == '2n3': # TODO: Untested
_3f = 11; _5f1 = 8; _5f2 = 2 # noqa
elif self.orientation == '2n3r': # TODO: Untested
_3f = 11; _5f1 = 2; _5f2 = 8 # noqa
elif self.orientation == '2n5': # TODO: Untested
_3f = 11; _5f1 = 8; _5f2 = 2 # noqa
elif self.orientation == '2n5r': # TODO: Untested
_3f = 11; _5f1 = 2; _5f2 = 8 # noqa
_3fold = self._3foldAxis[_3f]
_5fold_1 = self.getVertices()[_5f1]
_5fold_2 = self.getVertices()[_5f2]
plane1 = _normalizeVector(np.cross(_5fold_1, _3fold)) # cyan
plane2 = _normalizeVector(np.cross(_3fold, _5fold_2)) # magenta
plane3 = _normalizeVector(np.cross(_5fold_2, _5fold_1)) # green
if DEBUG:
logger.info("v0 %s" % _3fold)
logger.info("v1 %s" % _5fold_1)
logger.info("v2 %s" % _5fold_2)
logger.info("plane1 %s" % plane1)
logger.info("plane2 %s" % plane2)
logger.info("plane3 %s" % plane3)
return [_3fold, _5fold_1, _5fold_2], [plane1, plane2, plane3]