#
# The high-throughput toolkit (httk)
# Copyright (C) 2012-2015 Rickard Armiento
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero 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 Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import sys
from httk.core.geometry import is_point_inside_cell
from httk.core.httkobject import HttkObject, httk_typed_init, httk_typed_property
from httk.core.vectors import FracVector, FracScalar
from httk.core.basic import is_sequence
from httk.atomistic.cellutils import *
from httk.atomistic.cellshape import CellShape
from httk.atomistic.spacegrouputils import crystal_system_from_hall, lattice_system_from_hall
#TODO: Make this inherit CellShape
[docs]class Cell(HttkObject):
"""
Represents a cell (e.g., a unitcell, but also possibly just the basis vectors of a non-periodic system)
(The ability to represent the cell for a non-periodic system is also the reason this class is not called Lattice.)
"""
@httk_typed_init({'basis': (FracVector, 3, 3), 'lattice_system': str, 'orientation': int})
def __init__(self, basis, lattice_system, orientation=1):
"""
Private constructor, as per httk coding guidelines. Use Cell.create instead.
"""
self.basis = basis
self.orientation = orientation
self.lattice_system = lattice_system
self.niggli_matrix, self.orientation = basis_to_niggli_and_orientation(basis)
self.det = basis.det()
self.inv = basis.inv()
self._volume = abs(self.det)
self.metric = niggli_to_metric(self.niggli_matrix)
self.lengths, self.angles = niggli_to_lengths_and_angles(self.niggli_matrix)
self.lengths = [FracVector.use(x).simplify() for x in self.lengths]
self.angles = [FracVector.use(x).simplify() for x in self.angles]
_dummy, self.cosangles, self.sinangles = niggli_to_lengths_and_trigangles(self.niggli_matrix)
self.a, self.b, self.c = self.lengths
self.alpha, self.beta, self.gamma = self.angles
self.cosalpha, self.cosbeta, self.cosgamma = self.cosangles
self.sinalpha, self.sinbeta, self.singamma = self.sinangles
#def create(cls, basis=None, a=None, b=None, c=None, alpha=None, beta=None, gamma=None, volume=None, scale=None, niggli_matrix=None, orientation=1, lengths=None, angles=None, normalize=True):
[docs] @classmethod
def create(cls, cell=None, basis=None, metric=None, niggli_matrix=None,
a=None, b=None, c=None, alpha=None, beta=None, gamma=None,
lengths=None, angles=None, cosangles=None, scale=None,
scaling=None, volume=None, periodicity=None,
nonperiodic_vecs=None, orientation=1, hall=None,
lattice_system=None, eps=0):
"""
Create a new cell object,
cell: any one of the following:
- a 3x3 array with (in rows) the three basis vectors of the cell (a non-periodic system should conventionally use an identity matrix)
- a dict with a single key 'niggli_matrix' with a 3x2 array with the Niggli Matrix representation of the cell
- a dict with 6 keys, 'a', 'b', 'c', 'alpha', 'beta', 'gamma' giving the cell parameters as floats
scaling: free form input parsed for a scale.
positive value = multiply basis vectors by this value
negative value = rescale basis vectors so that cell volume becomes abs(value).
scale: set to non-None to multiply all cell vectors with this factor
volume: set to non-None if the basis vectors only give directions, and the volume of the cell should be this value (overrides scale)
periodicity: free form input parsed for periodicity
sequence: True/False for each basis vector being periodic
integer: number of non-periodic basis vectors
hall: giving the hall symbol makes it possible to determine the lattice system without numerical inaccuracy
lattice_system: any one of: 'cubic', 'hexagonal', 'tetragonal', 'orthorhombic', 'trigonal', 'triclinic', 'monoclinic', 'unknown'
"""
#print("Create cell:",cell,basis,angles, lengths,cosangles,a,b,c,alpha,beta,gamma)
if cell is not None:
return Cell.use(cell)
if basis is not None:
basis = FracVector.use(basis)
if angles is None and not (alpha is None or beta is None or gamma is None):
angles = [alpha, beta, gamma]
if lengths is None and not (a is None or b is None or c is None):
lengths = [a, b, c]
if cosangles is None and angles is not None:
cosangles = angles_to_cosangles(angles)
if basis is None and lengths is not None and cosangles is not None:
if hall is not None:
lattice_system = lattice_system_from_hall(hall)
else:
lattice_system = lattice_system_from_lengths_and_cosangles(lengths, cosangles)
basis = lengths_and_cosangles_to_conventional_basis(lengths, cosangles, lattice_system, orientation=orientation, eps=eps)
if niggli_matrix is not None:
niggli_matrix = FracVector.use(niggli_matrix)
if niggli_matrix is None and basis is not None:
niggli_matrix, orientation = basis_to_niggli_and_orientation(basis)
#if niggli_matrix is None and lengths is not None and cosangles is not None:
# niggli_matrix = lengths_and_cosangles_to_niggli(lengths, cosangles)
# niggli_matrix = FracVector.use(niggli_matrix)
#if niggli_matrix is None and lengths is not None and angles is not None:
# niggli_matrix = lengths_and_angles_to_niggli(lengths, angles)
# niggli_matrix = FracVector.use(niggli_matrix)
#if niggli_matrix is None and not (a is None or b is None or c is None or alpha is None or beta is None or gamma is None):
# niggli_matrix = lengths_and_angles_to_niggli([a, b, c], [alpha, beta, gamma])
# niggli_matrix = FracVector.use(niggli_matrix)
if basis is None:
raise Exception("cell.create: Not enough information to specify a cell given.")
if scaling is None and scale is not None:
scaling = scale
if scaling is not None and volume is not None:
raise Exception("Cell.create: cannot specify both scaling and volume!")
if volume is not None:
scaling = vol_to_scale(basis, volume)
if scaling is not None:
scaling = FracVector.use(scaling)
basis = (basis*scaling).simplify()
# Determination of the lattice system can only be made approximately, so it is recommended
# that it is given to the constructor if possible
if (lattice_system is None or lattice_system == 'unknown') and hall is not None:
lattice_system = lattice_system_from_hall(hall)
if lattice_system is None or lattice_system == 'unknown':
lattice_system = lattice_system_from_niggli(niggli_matrix)
if basis is None:
basis = FracVector.use(niggli_to_conventional_basis(niggli_matrix, lattice_system, orientation=orientation))
return cls(basis, lattice_system, orientation)
[docs] @classmethod
def use(cls, other):
if isinstance(other, Cell):
return other
else:
try:
if len(other) == 3:
return cls.create(basis=other)
elif len(other) == 2:
return cls.create(niggli_matrix=other)
elif len(other) == 1:
return cls.create(a=other[0], b=other[1], c=other[2], alpha=other[3], beta=other[4], gamma=other[5])
except Exception:
pass
raise Exception("Cell.use: do not know how to use an object of class:"+str(other.__class__))
@httk_typed_property(FracScalar)
def volume(self):
return self._volume
#@httk_typed_property((FracVector, 3, 3))
#def basis(self):
# return self._basis
#@httk_typed_property((FracVector, 3, 2))
#def niggli_matrix(self):
# return self._niggli_matrix
[docs] def scaling(self):
return -self.volume
[docs] def coords_reduced_to_cartesian(self, coords):
coords = FracVector.use(coords)
return coords*self.basis
[docs] def coordgroups_reduced_to_cartesian(self, coordgroups):
newcoordgroups = []
for coordgroup in coordgroups:
newcoordgroups += [self.coords_reduced_to_cartesian(coordgroup)]
return FracVector.stack(newcoordgroups)
[docs] def coords_cartesian_to_reduced(self, coords):
coords = FracVector.use(coords)
return coords*self.inv
[docs] def coordgroups_cartesian_to_reduced(self, coordgroups):
newcoordgroups = []
for coordgroup in coordgroups:
newcoordgroups += [self.coords_cartesian_to_reduced(coordgroup)]
return FracVector.stack(newcoordgroups)
[docs] def is_point_inside(self, cartesian_coord):
is_point_inside_cell(self.basis, cartesian_coord)
# # Type translators where convenient types and primitive types differ
# @classmethod
# def types_from_primitive(cls,data):
# if 'maxnorm_basis_g' in data:
# data['maxnorm_basis'] = FracVector(data['maxnorm_basis_g'],1000000000)
# del(data['maxnorm_basis_g'])
# if 'volume_g' in data:
# data['volume'] = FracVector(data['volume_g'],1000000000)
# del(data['volume_g'])
# return data
#
# @classmethod
# def types_to_primitive(cls,data):
# if 'maxnorm_basis' in data:
# adjustedmatrix = (data['maxnorm_basis_g']*1000000000).to_fractions()
# data['maxnorm_basis_g'] = [[int(x) for x in y] for y in adjustedmatrix]
# del(data['maxnorm_basis'])
# if 'volume_g' in data:
# data['volume'] = FracVector(data['volume_g'],1000000000)
# del(data['volume'])
# return data
[docs] def get_normalized_longestvec(self):
# Somewhat unusual normalization where the largest one element
# in the cell = 1, this way we avoid floating point operations for prototypes created from cell vectors
# (prototypes created from lengths and angles is another matter)
scale = self.normalization_scale.inv()
return Cell.create(basis=self.basis*scale.simplify())
[docs] def clean(self):
newbasis = self.basis.limit_denominator(5000000)
#new_niggli_matrix = self.niggli_matrix.limit_denominator(5000000)
return self.__class__(newbasis, lattice_system=self.lattice_system, orientation=self.orientation)
@property
def normalization_longestvec_scale(self):
"""
Get the factor with which a normalized version of this cell needs to be multiplied to reproduce this cell.
I.e. self = (normalization_scale)*self.get_normalized()
"""
# We use the somewhat unusual normalization where the largest one element
# in the cell = 1, this way we avoid floating point operations for prototypes created from cell vectors
# (prototypes created from lengths and angles is another matter)
#
c = self.basis
maxele = max(c[0, 0], c[0, 1], c[0, 2], c[1, 0], c[1, 1], c[1, 2], c[2, 0], c[2, 1], c[2, 2])
maxeleneg = max(-c[0, 0], -c[0, 1], -c[0, 2], -c[1, 0], -c[1, 1], -c[1, 2], -c[2, 0], -c[2, 1], -c[2, 2])
if maxeleneg > maxele:
scale = (-maxeleneg).simplify()
else:
scale = (maxele).simplify()
return scale
@property
def normalization_scale(self):
"""
"""
return FracScalar.create(pow(float(self.basis.det()), 0.33333333333333333333333))
[docs] def get_normalized(self):
scale = self.normalization_scale.inv()
return Cell.create(basis=self.basis*scale.simplify())
def __str__(self):
return "<Cell: %.8f * \n [" % (self.normalization_scale)+"\n ".join(["% .8f % .8f % .8f" % (x[0], x[1], x[2]) for x in self.get_normalized().basis.to_floats()])+"]>"
if __name__ == "__main__":
main()