#!/usr/bin/env python
##############################################################################
#
# diffpy.structure by DANSE Diffraction group
# Simon J. L. Billinge
# (c) 2006 trustees of the Michigan State University.
# All rights reserved.
#
# File coded by: Pavol Juhas
#
# See AUTHORS.txt for a list of people who contributed.
# See LICENSE_DANSE.txt for license information.
#
##############################################################################
"""Symmetry utility functions such as expansion of asymmetric unit, and
generation of positional constraints.
Attributes
----------
epsilon : float
Default tolerance for equality of 2 positions, also
used for identification of special positions.
stdUsymbols : list
Standard symbols denoting elements of anisotropic thermal
displacement tensor.
"""
import re
import sys
import numpy
from diffpy.structure.structureerrors import SymmetryError
from diffpy.utils._deprecator import build_deprecation_message, deprecated
base = "diffpy.structure"
removal_version = "4.0.0"
isSpaceGroupLatPar_deprecation_msg = build_deprecation_message(
base,
"isSpaceGroupLatPar",
"is_space_group_latt_parms",
removal_version,
)
isconstantFormula_deprecation_msg = build_deprecation_message(
base,
"isconstantFormula",
"is_constant_formula",
removal_version,
)
positionDifference_deprecation_msg = build_deprecation_message(
base,
"positionDifference",
"position_difference",
removal_version,
)
nearestSiteIndex_deprecation_msg = build_deprecation_message(
base,
"nearestSiteIndex",
"nearest_site_index",
removal_version,
)
equalPositions_deprecation_msg = build_deprecation_message(
base,
"equalPositions",
"equal_positions",
removal_version,
)
expandPosition_deprecation_msg = build_deprecation_message(
base,
"expandPosition",
"expand_position",
removal_version,
)
nullSpace_deprecation_msg = build_deprecation_message(
base,
"nullSpace",
"null_space",
removal_version,
)
# Constants ------------------------------------------------------------------
epsilon = 1.0e-5
stdUsymbols = ["U11", "U22", "U33", "U12", "U13", "U23"]
# ----------------------------------------------------------------------------
@deprecated(isSpaceGroupLatPar_deprecation_msg)
def isSpaceGroupLatPar(spacegroup, a, b, c, alpha, beta, gamma):
"""'diffpy.structure.isSpaceGroupLatPar' is deprecated and will be
removed in version 4.0.0.
Please use 'diffpy.structure.is_space_group_latt_parms' instead.
"""
return is_space_group_latt_parms(spacegroup, a, b, c, alpha, beta, gamma)
[docs]
def is_space_group_latt_parms(spacegroup, a, b, c, alpha, beta, gamma):
"""Check if space group allows passed lattice parameters.
Parameters
----------
spacegroup : SpaceGroup
Instance of `SpaceGroup`.
a, b, c, alpha, beta, gamma : float
`Lattice` parameters.
Return
------
bool
``True`` when lattice parameters are allowed by space group.
Note
----
Crystal system rules:
Benjamin, W. A., Introduction to crystallography, New York (1969), p.60.
"""
# crystal system rules
# ref: Benjamin, W. A., Introduction to crystallography,
# New York (1969), p.60
def check_triclinic():
return True
def check_monoclinic():
rv = (alpha == gamma == 90) or (alpha == beta == 90)
return rv
def check_orthorhombic():
return alpha == beta == gamma == 90
def check_tetragonal():
return a == b and alpha == beta == gamma == 90
def check_trigonal():
rv = (a == b == c and alpha == beta == gamma) or (a == b and alpha == beta == 90 and gamma == 120)
return rv
def check_hexagonal():
return a == b and alpha == beta == 90 and gamma == 120
def check_cubic():
return a == b == c and alpha == beta == gamma == 90
crystal_system_rules = {
"TRICLINIC": check_triclinic,
"MONOCLINIC": check_monoclinic,
"ORTHORHOMBIC": check_orthorhombic,
"TETRAGONAL": check_tetragonal,
"TRIGONAL": check_trigonal,
"HEXAGONAL": check_hexagonal,
"CUBIC": check_cubic,
}
rule = crystal_system_rules[spacegroup.crystal_system]
return rule()
# Constant regular expression used in isconstantFormula().
# isconstantFormula runs faster when regular expression is not
# compiled per every single call.
_rx_constant_formula = re.compile(r"[-+]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)??(/[-+]?\d+)?$")
@deprecated(isconstantFormula_deprecation_msg)
def isconstantFormula(s):
"""'diffpy.structure.isconstantFormula' is deprecated and will be
removed in version 4.0.0.
Please use 'diffpy.structure.is_constant_formula' instead.
"""
return is_constant_formula(s)
# Helper class intended for this module only:
class _Position2Tuple(object):
"""Create callable object that converts fractional coordinates to a
tuple of integers with given precision. For presision close to zero
it will return a tuples of double.
Note
----
Helper class intended for local use only.
Parameters
----------
eps : float
Cutoff for equivalent coordinates.
Attributes
----------
eps : float
Cutoff for equivalent coordinates. When two coordinates map to the
same tuple, they are closer than `eps`.
"""
def __init__(self, eps=None):
if eps is None:
eps = epsilon
# ensure self.eps has exact machine representation
self.eps = eps + 1.0
self.eps = self.eps - 1.0
# no conversions for very small eps
if self.eps == 0.0 or 1.0 / self.eps > sys.maxsize:
self.eps = 0.0
return
def __call__(self, xyz):
"""Convert array of fractional coordinates to a tuple.
Parameters
----------
xyz : Iterable
Fractional coordinates.
Return
------
tuple
Tuple of 3 float when `eps` is zero, otherwise tuple of 3 int.
"""
# no conversion case
if self.eps == 0.0:
tpl = tuple(xyz % 1.0)
return tpl
# here we convert to integer
tpl = tuple([int((xi - numpy.floor(xi)) / self.eps) for xi in xyz])
return tpl
# End of class _Position2Tuple
@deprecated(positionDifference_deprecation_msg)
def positionDifference(xyz0, xyz1):
"""'diffpy.structure.positionDifference' is deprecated and will be
removed in version 4.0.0.
Please use 'diffpy.structure.position_difference' instead.
"""
return position_difference(xyz0, xyz1)
[docs]
def position_difference(xyz0, xyz1):
"""Smallest difference between two coordinates in periodic lattice.
Parameters
----------
xyz0, xyz1 : array_like
Fractional coordinates.
Return
------
dxyz : numpy.ndarray
Smallest difference between two coordinates in periodic lattice
with ``0 <= dxyz <= 0.5``.
"""
dxyz = numpy.asarray(xyz0) - xyz1
# map differences to [0,0.5]
dxyz = dxyz - numpy.floor(dxyz)
mask = dxyz > 0.5
dxyz[mask] = 1.0 - dxyz[mask]
return dxyz
@deprecated(nearestSiteIndex_deprecation_msg)
def nearestSiteIndex(sites, xyz):
"""'diffpy.structure.nearestSiteIndex' is deprecated and will be
removed in version 4.0.0.
Please use 'diffpy.structure.nearest_site_index' instead.
"""
# we use box distance to be consistent with _Position2Tuple conversion
return nearest_site_index(sites, xyz)
[docs]
def nearest_site_index(sites, xyz):
"""Index of the nearest site to a specified position.
Parameters
----------
sites : array_like
List of coordinates.
xyz : array_like
Single position.
Return
------
int
Index of the nearest site.
"""
# we use box distance to be consistent with _Position2Tuple conversion
dbox = position_difference(sites, xyz).max(axis=1)
nearindex = numpy.argmin(dbox)
return nearindex
@deprecated(equalPositions_deprecation_msg)
def equalPositions(xyz0, xyz1, eps):
"""'diffpy.structure.equalPositions' is deprecated and will be
removed in version 4.0.0.
Please use 'diffpy.structure.equal_positions' instead.
"""
return equal_positions(xyz0, xyz1, eps)
[docs]
def equal_positions(xyz0, xyz1, eps):
"""Equality of two coordinates with optional tolerance.
Parameters
----------
xyz0, xyz1 : array_like
Fractional coordinates.
eps : float
Tolerance for equality of coordinates.
Return
------
bool
``True`` when two coordinates are closer than `eps`.
"""
# we use box distance to be consistent with _Position2Tuple conversion
dxyz = position_difference(xyz0, xyz1)
return numpy.all(dxyz <= eps)
@deprecated(expandPosition_deprecation_msg)
def expandPosition(spacegroup, xyz, sgoffset=[0, 0, 0], eps=None):
"""'diffpy.structure.expandPosition' is deprecated and will be
removed in version 4.0.0.
Please use 'diffpy.structure.expand_position' instead.
"""
return expand_position(spacegroup, xyz, sgoffset, eps)
[docs]
def expand_position(spacegroup, xyz, sgoffset=[0, 0, 0], eps=None):
"""Obtain unique equivalent positions and corresponding operations.
Parameters
----------
spacegroup : SpaceGroup
Instance of SpaceGroup.
xyz : list
Position to be expanded.
sgoffset : list, Optional
Offset of space group origin ``[0, 0, 0]``. Default is ``[0, 0, 0]``.
eps : float, Optional
Cutoff for equal positions, default is ``1.0e-5``.
Return
------
tuple
A tuple with ``(list of unique equivalent positions, nested
list of `SpaceGroups.SymOp` instances, site multiplicity)``.
"""
sgoffset = numpy.asarray(sgoffset, dtype=float)
if eps is None:
eps = epsilon
pos2tuple = _Position2Tuple(eps)
positions = []
site_symops = {} # position tuples with [related symops]
for symop in spacegroup.iter_symops():
# operate on coordinates in non-shifted spacegroup
pos = symop(xyz + sgoffset) - sgoffset
mask = numpy.logical_or(pos < 0.0, pos >= 1.0)
pos[mask] -= numpy.floor(pos[mask])
tpl = pos2tuple(pos)
if tpl not in site_symops:
pos_is_new = True
site_symops[tpl] = []
# double check if there is any position nearby
if positions:
nearpos = positions[nearest_site_index(positions, pos)]
# is it an equivalent position?
if equal_positions(nearpos, pos, eps):
# tpl should map to the same list as nearpos
site_symops[tpl] = site_symops[pos2tuple(nearpos)]
pos_is_new = False
if pos_is_new:
positions.append(pos)
# here tpl is inside site_symops
site_symops[tpl].append(symop)
# pos_symops is nested list of symops associated with each position
pos_symops = [site_symops[pos2tuple(p)] for p in positions]
multiplicity = len(positions)
return positions, pos_symops, multiplicity
@deprecated(nullSpace_deprecation_msg)
def nullSpace(A):
"""'diffpy.structure.nullSpace' is deprecated and will be removed in
version 4.0.0.
Please use 'diffpy.structure.null_space' instead.
"""
return null_space(A)
[docs]
def null_space(A):
"""Null space of matrix A."""
from numpy import linalg
u, s, v = linalg.svd(A)
# s may have smaller dimension than v
vnrows = numpy.shape(v)[0]
mask = numpy.ones(vnrows, dtype=bool)
mask[s > epsilon] = False
null_space = numpy.compress(mask, v, axis=0)
return null_space
def _find_invariants(symops):
"""Find a list of symmetry operations which contains identity.
Parameters
----------
symops : list of SymOp
Nested list of `SymOp` instances.
Return
------
list
List-item in symops which contains identity.
Raise
-----
ValueError
When identity is not found.
"""
invrnts = None
R0 = numpy.identity(3, dtype=float)
t0 = numpy.zeros(3, dtype=float)
for ops in symops:
for op in ops:
if numpy.all(op.R == R0) and numpy.all(op.t == t0):
invrnts = ops
break
if invrnts:
break
if invrnts is None:
emsg = "Could not find identity operation."
raise ValueError(emsg)
return invrnts
# ----------------------------------------------------------------------------
generator_site = "diffpy.symmetryutilities.GeneratorSite"
signedRatStr_deprecation_msg = build_deprecation_message(
generator_site,
"signedRatStr",
"convert_fp_num_to_signed_rational",
removal_version,
)
positionFormula_deprecation_msg = build_deprecation_message(
generator_site,
"positionFormula",
"position_formula",
removal_version,
)
UFormula_deprecation_msg = build_deprecation_message(
generator_site,
"UFormula",
"u_formula",
removal_version,
)
eqIndex_deprecation_msg = build_deprecation_message(
generator_site,
"eqIndex",
"eq_index",
removal_version,
)
[docs]
class GeneratorSite(object):
"""Storage of data related to a generator positions.
Parameters
----------
spacegroup : SpaceGroup
Instance of `SpaceGroup`.
xyz : array_like
Generating site. When `xyz` is close to special
position `self.xyz` will be adjusted.
Uij : array_like, Optional
Thermal factors at generator site. Yields `self.Uij`
after adjusting to spacegroup symmetry. Default is zeros.
sgoffset : list, Optional
Offset of space group origin ``[0, 0, 0]``. Default is ``[0, 0, 0]``.
eps : float, Optional
Cutoff for equal positions. Default is ``1.0e-5``.
Attributes
----------
xyz : numpy.ndarray
Fractional coordinates of generator site.
Uij : numpy.ndarray
Anisotropic thermal displacement at generator site.
sgoffset : numpy.ndarray
Offset of space group origin ``[0, 0, 0]``.
eps : float
Cutoff for equal positions.
eqxyz : list
List of equivalent positions.
eqUij : list
List of displacement matrices at equivalent positions.
symops : list
Nested list of operations per each `eqxyz`.
multiplicity : int
Generator site multiplicity.
Uisotropy : bool
Bool flag for isotropic thermal factors.
invariants : list
List of invariant operations for generator site.
null_space : numpy.ndarray
Null space of all possible differences of invariant
rotational matrices, this is a base of symmetry
allowed shifts.
Uspace : numpy.ndarray
3D array of independent components of U matrices.
pparameters : list
List of ``(xyz symbol, value)`` pairs.
Uparameters : list
List of ``(U symbol, value)`` pairs.
"""
Ucomponents = numpy.array(
[
[[1, 0, 0], [0, 0, 0], [0, 0, 0]],
[[0, 0, 0], [0, 1, 0], [0, 0, 0]],
[[0, 0, 0], [0, 0, 0], [0, 0, 1]],
[[0, 1, 0], [1, 0, 0], [0, 0, 0]],
[[0, 0, 1], [0, 0, 0], [1, 0, 0]],
[[0, 0, 0], [0, 0, 1], [0, 1, 0]],
],
dtype=float,
)
"""numpy.ndarray: 6x3x3 array of independent components of U
matrices."""
idx2Usymbol = {
0: "U11",
1: "U12",
2: "U13",
3: "U12",
4: "U22",
5: "U23",
6: "U13",
7: "U23",
8: "U33",
}
"""dict: Mapping of index to standard U symbol."""
def __init__(
self,
spacegroup,
xyz,
Uij=numpy.zeros((3, 3)),
sgoffset=[0, 0, 0],
eps=None,
):
if eps is None:
eps = epsilon
# just declare the members
self.xyz = numpy.array(xyz, dtype=float)
self.Uij = numpy.array(Uij, dtype=float)
self.sgoffset = numpy.array(sgoffset, dtype=float)
self.eps = eps
self.eqxyz = []
self.eqUij = []
self.symops = None
self.multiplicity = None
self.Uisotropy = False
self.invariants = []
self.null_space = None
self.Uspace = None
self.pparameters = []
self.Uparameters = []
# fill in the values
sites, ops, mult = expandPosition(spacegroup, xyz, sgoffset, eps)
invariants = _find_invariants(ops)
# shift self.xyz exactly to the special position
if mult > 1:
xyzdups = numpy.array([op(xyz + self.sgoffset) - self.sgoffset for op in invariants])
dxyz = xyzdups - xyz
dxyz = numpy.mean(dxyz - dxyz.round(), axis=0)
# recalculate if needed
if numpy.any(dxyz != 0.0):
self.xyz = xyz + dxyz
self.xyz[numpy.fabs(self.xyz) < self.eps] = 0.0
sites, ops, mult = expandPosition(spacegroup, self.xyz, self.sgoffset, eps)
invariants = _find_invariants(ops)
# self.xyz, sites, ops are all adjusted here
self.eqxyz = sites
self.symops = ops
self.multiplicity = mult
self.invariants = invariants
self._find_null_space()
self._find_pos_parameters()
self._find_u_space()
self._find_u_parameters()
self._find_eq_uij()
return
[docs]
def convert_fp_num_to_signed_rational(self, x):
"""Convert floating point number to signed rational
representation.
Possible fractional are multiples of 1/3, 1/6, 1/7, 1/9, if these
are not close, return `%+g` format.
Parameters
----------
x : float
Floating point number.
Return
------
str
Signed rational representation of `x`.
"""
s = "{:.8g}".format(x)
if len(s) < 6:
return "%+g" % x
den = numpy.array([3.0, 6.0, 7.0, 9.0])
nom = x * den
idx = numpy.where(numpy.fabs(nom - nom.round()) < self.eps)[0]
if idx.size == 0:
return "%+g" % x
# here we have fraction
return "%+.0f/%.0f" % (nom[idx[0]], den[idx[0]])
@deprecated(signedRatStr_deprecation_msg)
def signedRatStr(self, x):
"""'diffpy.structure.GeneratorSite.signedRatStr' is deprecated
and will be removed in version 4.0.0.
Please use 'diffpy.structure.GeneratorSite.convert_fp_num_to_signed_rational' instead.
"""
return self.convert_fp_num_to_signed_rational(x)
def _find_null_space(self):
"""Calculate `self.null_space` from `self.invariants`.
Try to represent `self.null_space` using small integers.
"""
R0 = self.invariants[0].R
Rdiff = [(symop.R - R0) for symop in self.invariants]
Rdiff = numpy.concatenate(Rdiff, axis=0)
self.null_space = null_space(Rdiff)
if self.null_space.size == 0:
return
# reverse sort rows of null_space rows by absolute value
key = tuple(numpy.fabs(numpy.transpose(self.null_space))[::-1])
order = numpy.lexsort(key)
self.null_space = self.null_space[order[::-1]]
# rationalize by the smallest element larger than cutoff
cutoff = 1.0 / 32
for row in self.null_space:
abrow = numpy.abs(row)
sgrow = numpy.sign(row)
# equalize items with round-off-equal absolute value
ii = abrow.argsort()
delta = 1e-8 * abrow[ii[-1]]
for k in ii[1:]:
if abrow[k] - abrow[k - 1] < delta:
abrow[k] = abrow[k - 1]
# find the smallest nonzero absolute element
jnz = numpy.flatnonzero(abrow > cutoff)
idx = jnz[abrow[jnz].argmin()]
row[:] = (sgrow * abrow) / sgrow[idx] / abrow[idx]
return
def _find_pos_parameters(self):
"""Find pparameters and their values for expressing
`self.xyz`."""
usedsymbol = {}
# parameter values depend on offset of self.xyz
txyz = self.xyz
# define txyz such that most of its elements are zero
for nvec in self.null_space:
idx = numpy.where(numpy.fabs(nvec) >= epsilon)[0][0]
varvalue = txyz[idx] / nvec[idx]
txyz = txyz - varvalue * nvec
# determine standard parameter name
vname = [s for s in "xyz"[idx:] if s not in usedsymbol][0]
self.pparameters.append((vname, varvalue))
usedsymbol[vname] = True
return
def _find_u_space(self):
"""Find independent U components with respect to invariant
rotations."""
n = len(self.invariants)
R6zall = numpy.tile(-numpy.identity(6, dtype=float), (n, 1))
R6zall_iter = numpy.split(R6zall, n, axis=0)
i6kl = (
(0, (0, 0)),
(1, (1, 1)),
(2, (2, 2)),
(3, (0, 1)),
(4, (0, 2)),
(5, (1, 2)),
)
for op, R6z in zip(self.invariants, R6zall_iter):
R = op.R
for j, Ucj in enumerate(self.Ucomponents):
Ucj2 = numpy.dot(R, numpy.dot(Ucj, R.T))
for i, kl in i6kl:
R6z[i, j] += Ucj2[kl]
Usp6 = null_space(R6zall)
# normalize Usp6 by its maximum component
mxcols = numpy.argmax(numpy.fabs(Usp6), axis=1)
mxrows = numpy.arange(len(mxcols))
Usp6 /= Usp6[mxrows, mxcols].reshape(-1, 1)
Usp6 = numpy.around(Usp6, 2)
# normalize again after rounding to get correct signs
mxcols = numpy.argmax(numpy.fabs(Usp6), axis=1)
Usp6 /= Usp6[mxrows, mxcols].reshape(-1, 1)
self.Uspace = numpy.tensordot(Usp6, self.Ucomponents, axes=(1, 0))
self.Uisotropy = len(self.Uspace) == 1
return
def _find_u_parameters(self):
"""Find Uparameters and their values for expressing
`self.Uij`."""
# permute indices as 00 11 22 01 02 12 10 20 21
diagorder = numpy.array((0, 4, 8, 1, 2, 5, 3, 6, 7))
Uijflat = self.Uij.flatten()
for Usp in self.Uspace:
Uspflat = Usp.flatten()
Uspnorm2 = numpy.dot(Uspflat, Uspflat)
permidx = next(i for i, x in enumerate(Uspflat[diagorder]) if x == 1)
idx = diagorder[permidx]
vname = self.idx2Usymbol[idx]
varvalue = numpy.dot(Uijflat, Uspflat) / Uspnorm2
self.Uparameters.append((vname, varvalue))
return
def _find_eq_uij(self):
"""Adjust `self.Uij` and `self.eqUij` to be consistent with
spacegroup."""
self.Uij = numpy.zeros((3, 3), dtype=float)
for i in range(len(self.Uparameters)):
Usp = self.Uspace[i]
varvalue = self.Uparameters[i][1]
self.Uij += varvalue * Usp
# now determine eqUij
for ops in self.symops:
# take first rotation matrix
R = ops[0].R
Rt = R.transpose()
self.eqUij.append(numpy.dot(R, numpy.dot(self.Uij, Rt)))
return
@deprecated(positionFormula_deprecation_msg)
def positionFormula(self, pos, xyzsymbols=("x", "y", "z")):
"""'diffpy.structure.GeneratorSite.positionFormula' is
deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.GeneratorSite.position_formula'
instead.
"""
return self.position_formula(pos, xyzsymbols)
@deprecated(UFormula_deprecation_msg)
def UFormula(self, pos, Usymbols=stdUsymbols):
"""'diffpy.structure.GeneratorSite.UFormula' is deprecated and
will be removed in version 4.0.0.
Please use 'diffpy.structure.GeneratorSite.u_formula' instead.
"""
return self.u_formula(pos, Usymbols)
@deprecated(eqIndex_deprecation_msg)
def eqIndex(self, pos):
"""'diffpy.structure.GeneratorSite.eqIndex' is deprecated and
will be removed in version 4.0.0.
Please use 'diffpy.structure.GeneratorSite.eq_index' instead.
"""
return self.eq_index(pos)
[docs]
def eq_index(self, pos):
"""Index of the nearest generator equivalent site.
Parameters
----------
pos : array_like
Fractional coordinates.
Return
------
int
Index of the nearest generator equivalent site.
"""
return nearest_site_index(self.eqxyz, pos)
# End of class GeneratorSite
# ----------------------------------------------------------------------------
[docs]
class ExpandAsymmetricUnit(object):
"""Expand asymmetric unit and anisotropic thermal displacement.
Parameters
----------
spacegroup : SpaceGroup
Instance of `SpaceGroup`.
corepos : array_like
List of positions in asymmetric unit,
it may contain duplicates.
coreUijs : numpy.ndarray, Optional
Thermal factors for `corepos`.
sgoffset : list, Optional
Offset of space group origin ``[0, 0, 0]``. Default is ``[0, 0, 0]``.
eps : float, Optional
Cutoff for duplicate positions. Default is ``1.0e-5``.
Attributes
----------
spacegroup : SpaceGroup
Instance of `SpaceGroup`.
corepos : array_like
List of positions in asymmetric unit,
it may contain duplicates.
coreUijs : numpy.ndarray
Thermal factors for `corepos`. Defaults to zeros.
sgoffset : numpy.ndarray
Offset of space group origin ``[0, 0, 0]``. Default to zeros.
eps : float
Cutoff for equivalent positions. Default is ``1.0e-5``.
multiplicity : list
Multiplicity of each site in `corepos`.
Uisotropy : list
Bool flags for isotropic sites in `corepos`.
expandedpos : list
List of equivalent positions per each site in `corepos`.
expandedUijs : list
List of thermal factors per each site in `corepos`.
"""
# By design Atom instances are not accepted as arguments to keep
# number of required imports low.
def __init__(self, spacegroup, corepos, coreUijs=None, sgoffset=[0, 0, 0], eps=None):
if eps is None:
eps = epsilon
# declare data members
self.spacegroup = spacegroup
self.corepos = corepos
self.coreUijs = None
self.sgoffset = numpy.array(sgoffset)
self.eps = eps
self.multiplicity = []
self.Uisotropy = []
self.expandedpos = []
self.expandedUijs = []
# obtain their values
corelen = len(self.corepos)
if coreUijs:
self.coreUijs = coreUijs
else:
self.coreUijs = numpy.zeros((corelen, 3, 3), dtype=float)
for cpos, cUij in zip(self.corepos, self.coreUijs):
gen = GeneratorSite(self.spacegroup, cpos, cUij, self.sgoffset, self.eps)
self.multiplicity.append(gen.multiplicity)
self.Uisotropy.append(gen.Uisotropy)
self.expandedpos.append(gen.eqxyz)
self.expandedUijs.append(gen.eqUij)
return
# End of class ExpandAsymmetricUnit
# Helper function for SymmetryConstraints class. It may be useful
# elsewhere therefore its name does not start with underscore.
pruneFormulaDictionary_deprecation_msg = build_deprecation_message(
base,
"pruneFormulaDictionary",
"prune_formula_dictionary",
removal_version,
)
@deprecated(pruneFormulaDictionary_deprecation_msg)
def pruneFormulaDictionary(eqdict):
"""'diffpy.structure.pruneFormulaDictionary' is deprecated and will
be removed in version 4.0.0.
Please use 'diffpy.structure.prune_formula_dictionary' instead.
"""
pruned = {}
for smb, eq in eqdict.items():
if not is_constant_formula(eq):
pruned[smb] = eq
return pruned
symmetry_constraints = "diffpy.symmetryutilities.SymmetryConstraints"
posparSymbols_deprecation_msg = build_deprecation_message(
symmetry_constraints,
"posparSymbols",
"pospar_symbols",
removal_version,
)
posparValues_deprecation_msg = build_deprecation_message(
symmetry_constraints,
"posparValues",
"pospar_values",
removal_version,
)
UparSymbols_deprecation_msg = build_deprecation_message(
symmetry_constraints,
"UparSymbols",
"upar_symbols",
removal_version,
)
UparValues_deprecation_msg = build_deprecation_message(
symmetry_constraints,
"UparValues",
"upar_values",
removal_version,
)
UFormulas_deprecation_msg = build_deprecation_message(
symmetry_constraints,
"UFormulas",
"u_formulas",
removal_version,
)
positionFormulas_deprecation_msg = build_deprecation_message(
symmetry_constraints,
"positionFormulas",
"position_formulas",
removal_version,
)
positionFormulasPruned_deprecation_msg = build_deprecation_message(
symmetry_constraints,
"positionFormulasPruned",
"position_formulas_pruned",
removal_version,
)
UFormulasPruned_deprecation_msg = build_deprecation_message(
symmetry_constraints,
"UFormulasPruned",
"u_formulas_pruned",
removal_version,
)
[docs]
class SymmetryConstraints(object):
"""Generate symmetry constraints for specified positions.
Parameters
----------
spacegroup : SpaceGroup
Instance of `SpaceGroup`.
positions : array_like
List of all positions to be constrained.
Uijs : array_like, Optional
List of U matrices for all constrained positions.
sgoffset : list, Optional
Offset of space group origin ``[0, 0, 0]``. Default is ``[0, 0, 0]``.
eps : float, Optional
Cutoff for duplicate positions. Default is ``1.0e-5``.
Attributes
----------
spacegroup : SpaceGroup
Instance of `SpaceGroup`.
positions : numpy.ndarray
All positions to be constrained.
Uijs : numpy.ndarray
Thermal factors for all positions. Defaults to zeros.
sgoffset : numpy.ndarray
Optional offset of space group origin ``[0, 0, 0]``.
eps : float
Cutoff for equivalent positions. Default is ``1.0e-5``.
corepos : list
List of of positions in the asymmetric unit.
coremap : dict
Dictionary mapping indices of asymmetric core positions
to indices of all symmetry related positions.
poseqns : list
List of coordinate formula dictionaries per each site.
Formula dictionary keys are from ``("x", "y", "z")`` and
the values are formatted as ``[[-]{x|y|z}%i] [{+|-}%g]``,
for example: ``x0``, ``-x3``, ``z7 +0.5``, ``0.25``.
pospars : list
List of ``(xyz symbol, value)`` pairs.
Ueqns : list
List of anisotropic atomic displacement formula
dictionaries per each position. Formula dictionary
keys are from ``('U11','U22','U33','U12','U13','U23')``
and the values are formatted as ``{[%g*][Uij%i]|0}``,
for example: ``U110``, ``0.5*U2213``, ``0``.
Upars : list
List of ``(U symbol, value)`` pairs.
Uisotropy : list
List of bool flags for isotropic thermal displacements.
"""
def __init__(self, spacegroup, positions, Uijs=None, sgoffset=[0, 0, 0], eps=None):
if eps is None:
eps = epsilon
# fill in data members
self.spacegroup = spacegroup
self.positions = None
self.Uijs = None
self.sgoffset = numpy.array(sgoffset)
self.eps = eps
self.corepos = []
self.coremap = {}
self.poseqns = None
self.pospars = []
self.Ueqns = None
self.Upars = []
self.Uisotropy = None
# handle list of lists returned by ExpandAsymmetricUnit
if len(positions) and isinstance(positions[0], list):
# concatenate lists before converting to Nx3 array
flatpos = sum(positions, [])
flatpos = numpy.array(flatpos, dtype=float).flatten()
self.positions = flatpos.reshape((-1, 3))
# otherwise convert to array
else:
flatpos = numpy.array(positions, dtype=float).flatten()
self.positions = flatpos.reshape((-1, 3))
# here self.positions should be a 2D numpy array
numpos = len(self.positions)
# adjust Uijs if not specified
if Uijs is not None:
self.Uijs = numpy.array(Uijs, dtype=float)
else:
self.Uijs = numpy.zeros((numpos, 3, 3), dtype=float)
self.poseqns = numpos * [None]
self.Ueqns = numpos * [None]
self.Uisotropy = numpos * [False]
# all members should be initialized here
self._find_constraints()
return
def _find_constraints(self):
"""Find constraints for positions and anisotropic displacements
`Uij`."""
numpos = len(self.positions)
# canonical xyzsymbols and Usymbols
xyzsymbols = [smbl + str(i) for i in range(numpos) for smbl in "xyz"]
Usymbols = [smbl + str(i) for i in range(numpos) for smbl in stdUsymbols]
independent = set(range(numpos))
for genidx in range(numpos):
if genidx not in independent:
continue
# it is a generator
self.coremap[genidx] = []
genpos = self.positions[genidx]
genUij = self.Uijs[genidx]
gen = GeneratorSite(self.spacegroup, genpos, genUij, self.sgoffset, self.eps)
# append new pparameters if there are any
gxyzsymbols = xyzsymbols[3 * genidx : 3 * (genidx + 1)]
for k, v in gen.pparameters:
smbl = gxyzsymbols["xyz".index(k)]
self.pospars.append((smbl, v))
gUsymbols = Usymbols[6 * genidx : 6 * (genidx + 1)]
for k, v in gen.Uparameters:
smbl = gUsymbols[stdUsymbols.index(k)]
self.Upars.append((smbl, v))
# search for equivalents inside indies
indies = sorted(independent)
for indidx in indies:
indpos = self.positions[indidx]
formula = gen.position_formula(indpos, gxyzsymbols)
# formula is empty when indidx is independent
if not formula:
continue
# indidx is dependent here
independent.remove(indidx)
self.coremap[genidx].append(indidx)
self.poseqns[indidx] = formula
self.Ueqns[indidx] = gen.u_formula(indpos, gUsymbols)
# make sure positions and Uijs are consistent with spacegroup
eqidx = gen.eq_index(indpos)
dxyz = gen.eqxyz[eqidx] - indpos
self.positions[indidx] += dxyz - dxyz.round()
self.Uijs[indidx] = gen.eqUij[eqidx]
self.Uisotropy[indidx] = gen.Uisotropy
# all done here
coreidx = sorted(self.coremap.keys())
self.corepos = [self.positions[i] for i in coreidx]
return
@deprecated(posparSymbols_deprecation_msg)
def posparSymbols(self):
"""'diffpy.structure.SymmetryConstraints.posparSymbols' is
deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.SymmetryConstraints.pos_parm_symbols' instead.
"""
return self.pos_parm_symbols()
[docs]
def pos_parm_symbols(self):
"""Return list of standard position parameter symbols."""
return [n for n, v in self.pospars]
@deprecated(posparValues_deprecation_msg)
def posparValues(self):
"""'diffpy.structure.SymmetryConstraints.posparValues' is
deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.SymmetryConstraints.pos_parm_values' instead.
"""
return self.pos_parm_values()
[docs]
def pos_parm_values(self):
"""Return list of position parameters values."""
return [v for n, v in self.pospars]
@deprecated(posparValues_deprecation_msg)
def positionFormulas(self, xyzsymbols=None):
"""'diffpy.structure.SymmetryConstraints.positionFormulas' is
deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.SymmetryConstraints.position_formulas' instead.
"""
return self.position_formulas(xyzsymbols)
@deprecated(positionFormulasPruned_deprecation_msg)
def positionFormulasPruned(self, xyzsymbols=None):
"""'diffpy.structure.SymmetryConstraints.positionFormulasPruned'
is deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.SymmetryConstraints.position_formulas_pruned' instead.
"""
return self.position_formulas_pruned(xyzsymbols)
@deprecated(UparSymbols_deprecation_msg)
def UparSymbols(self):
"""'diffpy.structure.SymmetryConstraints.UparSymbols' is
deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.SymmetryConstraints.u_parm_symbols' instead.
"""
return self.u_parm_symbols()
[docs]
def u_parm_symbols(self):
"""Return list of standard atom displacement parameter
symbols."""
return [n for n, v in self.Upars]
@deprecated(UparValues_deprecation_msg)
def UparValues(self):
"""'diffpy.structure.SymmetryConstraints.UparValues' is
deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.SymmetryConstraints.u_parm_values'
instead.
"""
return [v for n, v in self.Upars]
[docs]
def u_parm_values(self):
"""Return list of atom displacement parameters values."""
return [v for n, v in self.Upars]
@deprecated(UFormula_deprecation_msg)
def UFormulas(self, Usymbols=None):
"""'diffpy.structure.SymmetryConstraints.UFormulas' is
deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.SymmetryConstraints.u_formulas'
instead.
"""
return self.u_formulas(Usymbols)
@deprecated(UFormulasPruned_deprecation_msg)
def UFormulasPruned(self, Usymbols=None):
"""'diffpy.structure.SymmetryConstraints.UFormulasPruned' is
deprecated and will be removed in version 4.0.0.
Please use 'diffpy.structure.SymmetryConstraints.u_formulas_pruned'
instead.
"""
return self.u_formulas_pruned(Usymbols)
# End of class SymmetryConstraints
# ----------------------------------------------------------------------------
# basic demonstration
if __name__ == "__main__":
from diffpy.structure.spacegroups import sg100
site = [0.125, 0.625, 0.13]
Uij = [[1, 2, 3], [2, 4, 5], [3, 5, 6]]
g = GeneratorSite(sg100, site, Uij=Uij)
fm100 = g.position_formula(site)
print("g = GeneratorSite(sg100, %r)" % site)
print("g.positionFormula(%r) = %s" % (site, fm100))
print("g.pparameters =", g.pparameters)
print("g.Uparameters =", g.Uparameters)
print("g.UFormula(%r) =" % site, g.u_formula(site))