#
# This code is Copyright (C) 2015 The Cambridge Crystallographic Data Centre
# (CCDC) of 12 Union Road, Cambridge CB2 1EZ, UK and a proprietary work of CCDC.
# This code may not be used, reproduced, translated, modified, disassembled or
# copied, except in accordance with a valid licence agreement with CCDC and may
# not be disclosed or redistributed in any form, either in whole or in part, to
# any third party. All copies of this code made in accordance with a valid
# licence agreement as referred to above must contain this copyright notice.
#
# No representations, warranties, or liabilities are expressed or implied in the
# supply of this code by CCDC, its servants or agents, except where such
# exclusion or limitation is prohibited, void or unenforceable under governing
# law.
#
'''
The main class of the :mod:`ccdc.crystal` module is :class:`ccdc.crystal.Crystal`.
:class:`ccdc.crystal.Crystal` contains attributes and functions that relate
to crystallography. An example of a crystallographic attribute is the
:attr:`ccdc.crystal.Crystal.cell_volume`.
.. code-block:: python
>>> from ccdc.io import CrystalReader
>>> csd_crystal_reader = CrystalReader('CSD')
>>> first_csd_crystal = csd_crystal_reader[0]
>>> round(first_csd_crystal.cell_volume, 3)
769.978
'''
##################################################################################
import math
import collections
import functools
import itertools
import operator
from ccdc.molecule import Atom, Molecule, _file_object_factory
from ccdc.utilities import nested_class, _detect_format
from ccdc.utilities import _private_importer
with _private_importer() as pi:
pi.import_ccdc_module('UtilitiesLib')
pi.import_ccdc_module('MathsLib')
pi.import_ccdc_module('ChemistryLib')
pi.import_ccdc_module('ChemicalAnalysisLib')
pi.import_ccdc_module('FileFormatsLib')
pi.import_ccdc_module('MotifSearchLib')
pi.import_ccdc_module('PackingSimilarityLib')
try:
with _private_importer() as pi:
pi.import_ccdc_module('PddLib')
PddLib.CrystalStructureIdentifier()
except Exception:
_have_pdd = False
else:
_have_pdd = True
def _if_pdd_is_licensed(f):
'''Decorator to enable content only if PDD is licenced.'''
if _have_pdd:
return f
return None
##################################################################################
CellLengths = collections.namedtuple('CellLengths', ['a', 'b', 'c'])
CellAngles = collections.namedtuple('CellAngles', ['alpha', 'beta', 'gamma'])
[docs]class Crystal(object):
'''Represents a crystal.'''
_voids_volume_telemetry = 0
[docs] @nested_class('Crystal')
class ReducedCell(object):
'''The reduced cell of a crystal.'''
def __init__(self, _cell):
'''Private: initialiser.'''
self._reduced_cell = ChemistryLib.ReducedCell(_cell)
@property
def cell_lengths(self):
'''The cell lengths of the reduced cell.'''
return CellLengths(self._reduced_cell.a(), self._reduced_cell.b(), self._reduced_cell.c())
@property
def cell_angles(self):
'''The cell angles of the reduced cell.'''
return CellAngles(
self._reduced_cell.alpha().degrees(),
self._reduced_cell.beta().degrees(),
self._reduced_cell.gamma().degrees()
)
@property
def volume(self):
'''The volume of the reduced cell.'''
return self._reduced_cell.volume()
[docs] @nested_class('Crystal')
class HBond(Molecule.HBond):
'''A crystallographic hydrogen bond.'''
def __init__(self, _view, _contact):
'''Private.'''
self._view = _view
self._contact = _contact
@property
def type(self):
'''The type of contact.'''
return Molecule.Contact._type_map[self._contact.type()]
@property
def intermolecular(self):
'''Whether the contact is inter- or intra-molecular.'''
return self._contact.intermolecular()
@property
def symmetry_operators(self):
'''The symmetry operators by which the atoms of the contact are related.'''
def _f(at):
base = self._view.base_asymmetric_unit_atom(at._atom)
op = ChemistryLib.atom_atom_symmetry_relation(
self._view.crystal_structure(), base, at._atom
)
if op:
return op.to_string()
return ''
return [_f(a) for a in self.atoms]
[docs] @nested_class('Crystal')
class MillerIndices(object):
'''Represents the family of planes subtended by Miller indices.'''
def __init__(self, h, k, l, crystal, _cell=None):
if h == k == l == 0:
raise TypeError('Miller indices may not be uniformly 0')
self._miller_indices = ChemistryLib.MillerIndices(h, k, l)
self.crystal = crystal
if crystal is not None:
self._cell = crystal._crystal.cell()
else:
self._cell = _cell
def __repr__(self):
'''String representation.'''
return 'MillerIndices(%d, %d, %d)' % (
self._miller_indices.h(), self._miller_indices.k(), self._miller_indices.l()
)
def __hash__(self):
return hash(self._miller_indices)
@property
def hkl(self):
'''The indices.'''
return (self._miller_indices.h(), self._miller_indices.k(), self._miller_indices.l())
@property
def plane(self):
'''The plane intersecting the unit cell.'''
from ccdc.descriptors import GeometricDescriptors
return GeometricDescriptors.Plane(0, 0, _plane=self._cell.hklplane(self._miller_indices))
@property
def order(self):
'''The order for improper Miller indices.'''
return self._miller_indices.order()
@property
def proper(self):
'''The proper, relatively prime Miller indices.'''
proper = self._miller_indices.proper()
return Crystal.MillerIndices(
proper.h(), proper.k(), proper.l(), self.crystal
)
@property
def d_spacing(self):
try:
hplane = self._cell.hklplane(self._miller_indices)
return hplane.distance()
except RuntimeError:
pass
[docs] @nested_class('Crystal')
class Disorder():
'''A class to represent disorder in the crystal structure.'''
[docs] @nested_class('Disorder')
class Group():
'''Represents a disorder group in a disorder assembly.'''
def __init__(self, csv, assembly, group):
'''Initialize the Group object'''
self._csv = csv
self._assembly_index = assembly
self._group_index = group
[docs] def activate(self):
'''Set this disorder group as active.'''
self._csv.set_disorder_assembly_group(
self._assembly_index,
self._group_index)
@property
def id(self):
'''The identifier of this disorder group.'''
return self._group_index
@property
def occupancy(self):
'''The occupancy of this disorder group.'''
return self._csv.disorder_group_occupancy(
self._assembly_index,
self._group_index)
@property
def atoms(self):
'''The list of atoms in this disorder group.
:returns: A list of :class:`ccdc.molecule.Atom` objects.
'''
return [
Atom(_atom=atom)
for atom in self._csv.disorder_group_atoms(
self._assembly_index, self._group_index)
]
[docs] @nested_class('Disorder')
class Assembly():
'''Represents a disorder assembly.'''
def __init__(self, csv, assembly):
'''Initialize the Assembly object'''
self._csv = csv
self._assembly_index = assembly
self._groups = tuple(
Crystal.Disorder.Group(csv, assembly, group)
for group in range(csv.disorder_assembly_ngroups(assembly)))
@property
def id(self):
'''The identifier of the disorder assembly.'''
return self._assembly_index
@property
def groups(self):
'''The disorder groups in this assembly.
:returns: A tuple of :class:`ccdc.crystal.Crystal.Disorder.Group` objects.
'''
return self._groups
@property
def active(self):
'''The active disorder group of this assembly.
:returns: A :class:`ccdc.crystal.Crystal.Disorder.Group` object.
'''
group_index = self._csv.disorder_assembly_current_group(self._assembly_index)
return self._groups[group_index]
def __init__(self, csv):
'''Initialize the :class:`ccdc.crystal.Crystal.Disorder` object'''
self._csv = csv
self._assemblies = tuple(
Crystal.Disorder.Assembly(csv, assembly)
for assembly in range(self._csv.ndisorder_assemblies()))
@property
def assemblies(self):
'''The disorder assemblies in the crystal structure.
:returns: A tuple of :class:`ccdc.crystal.Crystal.Disorder.Assembly` objects.
'''
return self._assemblies
@property
def is_suppressed(self):
'''Return True if disorder is suppressed, False if disorder is
analysed.
If the disorder is suppressed, all of the minority occupancy, that
is, suppressed, atoms are contained in a single disorder assembly.
(The majority occupancy atoms are always bonded and so are not
included in the disorder object.) The disorder assembly contains
two disorder groups. The first one includes all the suppressed
atoms, and the second none of them. The second group is by default
the active group, which describes a molecule with just the majority
occupancy atoms and none of the minority occupancy atoms.
If the disorder is not suppressed, it is fully represented.
Assemblies exist for every independent area of disorder in the
crystal. And multiple groups exist in each assembly.
'''
return self._csv.has_suppressed_disorder()
@property
def combinations(self):
'''Yield combination of disorder groups.
Note that the crystal is updated with the yielded disorder groups.
'''
# Save current active groups
current_groups = [assembly.active.id for assembly in self.assemblies]
# Iterate over all combinations
group_index_combinations = itertools.product(*[
list(range(len(assembly.groups)))
for assembly in self.assemblies])
try:
for group_index_combination in group_index_combinations:
active_groups = tuple(
self.assemblies[assembly].groups[group]
for assembly, group in enumerate(group_index_combination))
for active_group in active_groups:
active_group.activate()
yield active_groups
finally:
# Re-activate saved active groups
for assembly, group in enumerate(current_groups):
self.assemblies[assembly].groups[group].activate()
def __init__(self, crystal, identifier):
'''Create crystal using a ChemistryLib.CrystalStructure and an identifier.'''
self._crystal = crystal
self._identifier = identifier
self._csv = ChemistryLib.CrystalStructureView_instantiate(crystal)
# For crystal structure with suppressed disorder, make second disorder
# group the default as .molecule takes into account disorder selection
if self._csv.has_suppressed_disorder():
self._csv.set_disorder_assembly_group(0, 1)
def __eq__(self, other):
'''Equality of entries.'''
return isinstance(other, Crystal) and self.identifier == other.identifier
def __hash__(self):
'''So entries may be placed in dictionaries.'''
return hash(self.identifier)
@property
def identifier(self):
'''The string identifier of the crystal, e.g. 'ABEBUF'.'''
return self._identifier
@identifier.setter
def identifier(self, value):
self._identifier = value
@property
def reduced_cell(self):
'''The reduced cell of the crystal.'''
return Crystal.ReducedCell(self._crystal.cell())
@property
def orthogonal_to_fractional (self):
'''Return the transformation mapping orthogonal to fractional coordinates in the crystal cell.'''
o2f = self._crystal.cell().o2f()
transformation = MathsLib.Transformation(
o2f, MathsLib.vector_3d(0,0,0))
return Molecule.Transformation(_transformation=transformation)
@property
def fractional_to_orthogonal (self):
'''Return the transformation mapping fractional to orthogonal coordinates in the crystal cell.'''
f2o = self._crystal.cell().f2o()
transformation = MathsLib.Transformation(
f2o, MathsLib.vector_3d(0,0,0))
return Molecule.Transformation(_transformation=transformation)
@property
def cell_lengths(self):
'''Tuple containing the cell axis lengths: (a, b, c).
The returned value may be addressed by index or key:
>>> l = CellLengths(1, 2, 3)
>>> l[0]
1
>>> l.b
2
'''
return CellLengths(
self._crystal.cell().a(),
self._crystal.cell().b(),
self._crystal.cell().c()
)
#@property
#def uncertain_cell_lengths(self):
# '''Tuple containing :class:`ccdc.utilities.UncertainValue` for the lengths.'''
# try:
# return CellLengths(
# UncertainValue(self._crystal.cell().a(), self._crystal.cell().a_uncertainty()),
# UncertainValue(self._crystal.cell().b(), self._crystal.cell().b_uncertainty()),
# UncertainValue(self._crystal.cell().c(), self._crystal.cell().c_uncertainty())
# )
# except:
# pass
@property
def cell_angles(self):
'''Tuple containing the cell angles; (alpha, beta, gamma).
Note that the angles are reported in degrees.
The returned value may be addressed by index or by key:
>>> a = CellAngles(90.0, 45.0, 33.3)
>>> a.alpha
90.0
>>> a[1]
45.0
'''
return CellAngles(
self._crystal.cell().alpha().degrees(),
self._crystal.cell().beta().degrees(),
self._crystal.cell().gamma().degrees()
)
@property
def has_disorder(self):
'''Whether or not the crystal has disorder.'''
return self._crystal.ndisorder_assemblies() > 0
@property
def lattice_centring(self):
'''The lattice centring of this crystal.'''
return ChemistryLib.Spacegroup_centring_text().text(
self._crystal.cell().spacegroup().centring()
)
#@property
#def uncertain_cell_angles(self):
# '''Tuple containing :class:`ccdc.utilities.UncertainValue` for the cell angles.'''
# try:
# return CellAngles(
# UncertainValue(self._crystal.cell().alpha().degrees(), self._crystal.cell().alpha_uncertainty()),
# UncertainValue(self._crystal.cell().beta().degrees(), self._crystal.cell().beta_uncertainty()),
# UncertainValue(self._crystal.cell().gamma().degrees(), self._crystal.cell().gamma_uncertainty())
# )
# except:
# pass
@property
def z_value(self):
'''The number of molecules in the unit cell.'''
try:
return self._crystal_info.z_value() # pylint: disable=E1101
except AttributeError:
return None
@property
def z_prime(self):
'''The number of molecules in the asymmetric unit.'''
try:
return self._crystal_info.z_prime() # pylint: disable=E1101
except AttributeError:
return None
@property
def spacegroup_symbol(self):
'''The space group symbol of the crystal.'''
return self._crystal.cell().spacegroup().name()
@property
def crystal_system(self):
'''The space group system of the crystal.'''
return self._crystal.cell().spacegroup().crystal_system_text(
self._crystal.cell().spacegroup().crystal_system() )
@property
def spacegroup_number_and_setting(self):
'''The number in international tables and setting of the crystal's space group.
>>> from ccdc.io import CrystalReader
>>> csd_crystal_reader = CrystalReader('CSD')
>>> crystal = csd_crystal_reader.crystal('AABHTZ')
>>> crystal.spacegroup_number_and_setting
(2, 1)
Non standard spacegroup numbers, those above 230, will be returned with setting number 0.
Unrecognised spacegroups will raise a RuntimeError.
'''
num = self._crystal.cell().spacegroup().number()
if num > 230:
return (num, 0)
return ChemistryLib.InternationalTables.number_and_setting(
[self._crystal.cell().spacegroup().symmop(i) for i in range(self._crystal.cell().spacegroup().nsymmops())]
)
[docs] @staticmethod
def symmetry_operator_description(symmetry_operator):
'''A textual description of the symmetry operator.'''
s = symmetry_operator.replace(' ', '')
symm = ChemistryLib.CrystalSymmetryOperator(s)
return symm.detailed_description()
@property
def symmetry_operators(self):
'''The symmetry operators pertaining to this crystal.
:return: a tuple of string representations of the symmetry operators
'''
return tuple(
self._crystal.cell().spacegroup().symmop(i).to_string()
for i in range(self._crystal.cell().spacegroup().nsymmops())
)
[docs] @staticmethod
def symmetry_rotation(operator):
'''The rotational component of the symmetry operator.
:param operator: a string representation of a symmetry operator
:returns: a tuple of 9 integer values representing the 3x3 rotation matrix of the operator
'''
m = ChemistryLib.CrystalSymmetryOperator(operator.replace(' ', ''))
r = m.rotation()
return (
r.row(1).x(), r.row(1).y(), r.row(1).z(),
r.row(2).x(), r.row(2).y(), r.row(2).z(),
r.row(3).x(), r.row(3).y(), r.row(3).z()
)
[docs] @staticmethod
def symmetry_translation(operator):
'''The translational component of the symmetry operator.
:param operator: a string representation of a symmetry operator
:returns: a tuple of three floats representing the translational component of the operator
'''
m = ChemistryLib.CrystalSymmetryOperator(operator.replace(' ', ''))
t = m.trans()
return (t.x(), t.y(), t.z())
[docs] def atoms_on_special_positions(self, symmetry_operator=None):
'''The tuple of atoms lying on symmetry axes.
:param symmetry_operator: a symmetry operator, or ``None`` to get atoms on any symmetry axis.
:returns: a set of :class:`ccdc.molecule.Atom` instances.
'''
mol = self.molecule
atoms = mol.atoms
def _close_enough(a, b):
if a.coordinates is not None:
return all(
abs(a.coordinates[i] - b.coordinates[i]) < 1e-6 for i in range(3)
)
def _special_atoms(op):
symm = self.symmetric_molecule(op)
symm_ats = symm.atoms
unmoved = set(a for a, b in zip(atoms, symm_ats) if _close_enough(a, b))
return unmoved
if symmetry_operator is not None:
return _special_atoms(symmetry_operator)
else:
return functools.reduce(operator.__or__, tuple(_special_atoms(op) for op in self.symmetry_operators)[1:], set())
[docs] def symmetric_molecule(self, symmop, translate=None, force=False):
'''Generate a symmetry related copy of the molecule.
This method may be used to build multi-molecular crystals for visualisation and other purposes.
:param symmop: a string representation of the symmetry operation, such as '-x,1/2+y,1/2-z' representing a 2-fold screw axis.
:param translate: a sequence of three integers representing the translational component to be applied to the symmetry operation. If this is None, then the symmetry operator will contain the unit cell translation parameters, e.g. of the form '1-x,2-y,-1+z'.
:param force: a boolean value to allow symmetry operators not supported by the crystal to be applied.
:returns: a :class:`ccdc.molecule.Molecule` derived from this crystal with the symmetry operation applied.
:raises: TypeError if the symmetry operator is not in the symmetry operators of this crystal and force is not True.
'''
s = symmop.replace(' ', '')
if translate is None:
symm = ChemistryLib.CrystalSymmetryTransformation(s)
s = symm.symmop().to_string()
else:
t = MathsLib.int_vector_3d(int(translate[0]), int(translate[1]), int(translate[2]))
symm = ChemistryLib.CrystalSymmetryTransformation(
ChemistryLib.CrystalSymmetryOperator(s), t
)
if s not in self.symmetry_operators and not force:
raise TypeError('The symmetry operator %s is not applicable to this crystal' % symmop)
m = ChemistryLib.create_editable_molecule(
self.molecule._molecule, symm, self._crystal.cell()
)
m = m.create_editable_molecule()
m = Molecule(self.identifier, _molecule=m)
m._cell = self._crystal.cell()
return m
[docs] def centre_molecule(self):
'''Ensures that the molecule of the crystal lies within the unit cell.'''
self._crystal.move_molecules_to_cell_centre()
@property
def molecule(self):
'''The underlying molecule.
Note that a molecule can have several components, and, where present,
this takes into account the disorder selection.
'''
if self.has_disorder:
m = self._csv.super_molecule()
else:
m = self._crystal.editable_molecule()
m = ChemistryLib.Molecule3d_instantiate(m)
m = Molecule(self.identifier, _molecule=m).copy()
m.remove_atoms(a for a in m.atoms if a.label.endswith('?'))
m._cell = self._crystal.cell()
return m
@molecule.setter
def molecule(self, molecule):
'''Set the molecule.'''
edmol = molecule._molecule.create_editable_molecule()
for i in range(edmol.natoms()):
a = edmol.atom(i)
if a.site():
a.site().set_orth(a.site().orth(), self._crystal.cell())
self._crystal.set_editable_molecule(edmol)
@property
def disordered_molecule(self):
'''The underlying molecule with disordered atoms represented.
'''
if self.has_disorder and not self.disorder.is_suppressed:
# With a structure with analysed disorder, suppress the atoms so
# that there are no more bonds than there should be.
# This is needed as disorder is now analysed when CIF is read.
cs = self._crystal.clone()
ChemicalAnalysisLib.DisorderSuppresser().suppress_atoms(cs)
m = cs.editable_molecule()
else:
m = self._crystal.editable_molecule()
m = Molecule(self.identifier, _molecule=m).copy()
m._cell = self._crystal.cell()
return m
@property
def asymmetric_unit_molecule(self):
'''The molecular representation of the asymmetric unit.'''
csv = ChemistryLib.CrystalStructureView_instantiate(self._csv)
csv.asymmetric_unit()
mol = ChemistryLib.Molecule3d_instantiate()
for i in range(csv.nmolecules()):
mol.add(csv.molecule(i))
return Molecule(self.identifier, _molecule=mol)
[docs] def hbonds(self, intermolecular='Any', distance_range=(-5.0, 0.0), angle_tolerance=120.0, vdw_corrected=True, require_hydrogens=True, path_length_range=(4, 999), hbond_criterion=None, unique=True):
'''The collection of crystallographic hydrogen bonds in this crystal structure.
A crystallographic hydrogen bond may span two symmetry related molecules of the crystal.
The definition of a hydrogen bond (*D-H...A*) is a contact meeting the following criteria:
* The donor (*D*) must be a nitrogen, oxygen or sulphur atom covalently bound to at least one hydrogen.
* The acceptor (*A*) must be a nitrogen, oxygen, sulphur or halogen with at least one available lone
pair (*e.g.*, pyramidal trigonal nitrogen is regarded as an acceptor but planar trigonal nitrogen is not).
* The distance must be within the specified range (`distance_range`)
relative to the sum of van der Waals Radii of the atoms involed
in the H-bond (`vdw_corrected` is `True`) or within the absolute
specified range (`vdw_corrected` is `False`) in Angstroms. If
`require_hydrogens` is `True`, the distance constraint refers to
the *H...A* distance; otherwise it refers to the *D...A* distance.
* The *D-H...A* angle must be within the specified `angle_tolerance` range.
* The contact may be intermolecular and/or intramolecular involving donor and acceptor atoms separated with
covalent bonds within the specified path length.
If hydrogens are not required, they will not appear in the returned :attr:`ccdc.molecule.Molecule.HBond.atoms`.
:param intermolecular: 'Intramolecular', 'Intermolecular' or 'Any'.
:param distance_range: Minimum and maximum distances considered acceptable for a hydrogen bond to be formed (Angstroms).
:param vdw_corrected: Whether the distances are relative to the Van der Waals radius of the atoms.
:param angle_tolerance: Acceptable range for a hydrogen Donor-Hydrogen-Acceptor angle (degrees).
:param require_hydrogens: Whether hydrogens are necessary for the hydrogen bond.
:param path_length_range: Minimum and maximum values for length of path between donor and acceptor atom
:param hbond_criterion: an instance of :class:`ccdc.molecule.Molecule.HBondCriterion` or ``None``. If not ``None`` the definitions there will override any other specified argument.
:param unique: bool. If True, the crystallographically unique HBonds will be returned; if ``False`` the full set of HBonds will be returned. These may include some symmetry copies.
:returns: a tuple of :class:`ccdc.crystal.Crystal.HBond`
'''
if hbond_criterion is not None:
hbc = hbond_criterion._hbc.hbond_criterion()
else:
if intermolecular.lower().startswith('intra'):
inter = ChemistryLib.ContactCriterion.INTRAMOLECULAR
elif intermolecular.lower().startswith('inter'):
inter = ChemistryLib.ContactCriterion.INTERMOLECULAR
else:
inter = ChemistryLib.ContactCriterion.ANY
angle = MathsLib.Angle(angle_tolerance, MathsLib.Angle.DEGREES)
hbc = ChemistryLib.HBondCriterion(
max(distance_range), min(distance_range), angle,
vdw_corrected, not require_hydrogens,
inter,
min(path_length_range), max(path_length_range)
)
view = self._csv
ret = []
def _base_atom(at):
return view.base_asymmetric_unit_atom(at._atom)
for contact in view.find_base_contacts(hbc):
hb = Crystal.HBond(_view=view, _contact=contact)
for h in ret:
if (
_base_atom(hb.atoms[0]) == _base_atom(h.atoms[-1]) and
_base_atom(hb.atoms[-1]) == _base_atom(h.atoms[0]) and
(unique or hb.symmetry_operators[0] == hb.symmetry_operators[-1] == 'x,y,z')
):
break
else:
ret.append(hb)
return tuple(ret)
@property
def formula(self):
'''Return the chemical formula of the molecule in the crystal.'''
csv = ChemistryLib.CrystalStructureView_instantiate(self._crystal)
form = ChemistryLib.CrystalChemicalFormula(csv)
return str(form.sum_formula())
@property
def cell_volume(self):
'''Volume of the unit cell.'''
return self._crystal.cell().volume()
@property
def calculated_density(self):
'''Calculated density of the crystal.'''
if self.z_value is None or self.z_prime is None or self.z_prime == 0.0:
return None
try:
effective_z_prime = math.ceil(self.z_prime) if self.z_prime > 1 else 1
return ChemistryLib.CrystalDensityCalculator().calculate_experimental_density(
self.formula, self.z_value, self.cell_volume
) / effective_z_prime
except RuntimeError:
pass
@property
def packing_coefficient(self):
'''The packing coefficient of the crystal.
Measures the proportion of the unit cell occupied by atoms. It is a
fraction between zero and one; going from unoccupied to completely
filled.
'''
if self._crystal.cell().is_valid():
try:
return ChemicalAnalysisLib.CrystalStructureProperties().packing_coefficient(self._crystal)
except:
pass
# @property
# def asymmetric_unit_mass(self):
# '''The combined mass of the symmetry independent atoms of the crystal.
#
# This property ignores atoms that are do not have coordinates.
# '''
# # Make sure that the CrystalStructureProperties().unit_cell_mass is
# # really returning the asymmetric unit mass. Which it was doing at
# # 13/05/2014.
# return ChemicalAnalysisLib.CrystalStructureProperties().unit_cell_mass(self._crystal)
[docs] def void_volume(self, probe_radius=1.2, grid_spacing=0.7, mode='contact'):
'''Determine the void volume of the crystal.
>>> from ccdc.io import EntryReader
>>> entry_reader = EntryReader('CSD')
>>> abawop_crystal = entry_reader.crystal('ABAWOP')
>>> round(abawop_crystal.void_volume(), 2)
13.99
:param probe_radius: float, size of the probe
:param grid_spacing: float, fineness of the grid on which the calculation is made
:param mode: either 'accessible' or 'contact' according to whether the
centre of the probe or the whole probe must be accomodated.
:returns: void volume as a percentage of the unit cell volume
'''
if type(self)._voids_volume_telemetry == 0:
UtilitiesLib.ccdc_voids_volume_telemetry()
type(self)._voids_volume_telemetry = 1
vf = ChemicalAnalysisLib.CrystalVoidsFinder()
atoms = vf.calculate_atoms(self._csv)
if not atoms:
raise RuntimeError('No atoms in the crystal')
if mode.lower().startswith('access'):
mode = vf.ACCESSIBLE_SURFACE
else:
mode = vf.CONTACT_SURFACE
vol = vf.find_voids_volume(self._csv, probe_radius, grid_spacing, mode)
return vol
[docs] def assign_bonds(self):
'''Detect and assign bond types in the crystal.
This function will ignore any existing bond information and will use
the geometry of the crystal to detect and assign bond types.
'''
bd = ChemicalAnalysisLib.CrystalBondDetector()
bd.detect_bonds(self._crystal)
assigner = ChemicalAnalysisLib.MoleculeAllBondTypesAssigner()
assigner.assign(self._crystal.editable_molecule())
[docs] def remove_hydrogens(self):
'''Remove all hydrogen atoms from this crystal'''
ChemistryLib.remove_hydrogens(self._crystal.editable_molecule())
[docs] def add_hydrogens(self, mode='all', add_sites=True):
'''Add hydrogen atoms to the crystal.
This method adds hydrogens to a crystal (as opposed to a molecule).
Unlike the associated Molecule method :func:`~ccdc.Molecule.add_hydrogens` this version
adds to the crystal and is more crystallographically aware (trying to account for
possible disorder and also the crystalline environment when adding sites.)
:param mode: 'all' to generate all hydrogens (throws away existing
hydrogens) or 'missing' to generate hydrogens deemed to
be missing.
:param add_sites: 'True' to generate 3D coordinates for hydrogens (default) 'False' to only generate siteless atoms
:raises: RuntimeError if any heavy atom has no site.
:raises: RuntimeError if any atoms are of unknown type.
'''
if mode.lower() == 'all':
self.remove_hydrogens()
ChemistryLib.molecule_add_missing_hydrogens(self._crystal, add_sites)
[docs] def packing_shell(self, packing_shell_size=15):
'''Create a packing shell of the crystal.
This method is available only to CSD-Materials and CSD-Enterprise users.
:param packing_shell_size: the required number of molecules in the packing shell
:return: a :class:`ccdc.molecule.Molecule` containing packing_shell_size replicas of the crystal.
'''
calc = PackingSimilarityLib.MolecularPackingShellCalculator()
calc.set_shell_size(packing_shell_size)
view = ChemistryLib.CrystalStructureView.instantiate(self._crystal)
mols = calc.create_packing_shell(view)
shell = Molecule(self.identifier)
for m in mols:
shell.add_molecule(Molecule('', _molecule=m))
return shell
[docs] def molecular_shell(self, distance_type="actual",
distance_range=(0.0, 3.0),
atom_selection=[],
include_central=False):
"""
Return a molecule containing the atoms within a distance cutoff.
A subset of atoms to base the expansion on can be provided using the
atom_selection argument.
If the distance type is VdW then the min to max range is relative to
the sum of the vdW radii.
:param distance_type: 'vdw' or 'actual'.
:param distance_range: Minimum and maximum distances for the contact distance range.
:param atom_selection: list of atoms to base the expansion on.
:param include_central: if True the returned shell will contain the central molecules,
otherwise they are excluded.
:returns: :class:`ccdc.molecule.Molecule` containing the atoms within
the distance cutoff.
"""
crystal_view = ChemistryLib.CrystalStructureView.instantiate(self._crystal)
contact_criterion = ChemistryLib.CombinedCriterion()
contact_criterion.set_only_strongest(True)
contact_criterion.set_min_tolerance(min(distance_range))
contact_criterion.set_max_tolerance(max(distance_range))
if distance_type.lower() == "vdw":
contact_criterion.set_vdw_corrected(True)
else:
contact_criterion.set_vdw_corrected(False)
packing_molecule = ChemistryLib.Molecule3d_instantiate()
core_mols = list(crystal_view.molecule(i) for i in range(crystal_view.nmolecules()))
if include_central:
for mol in core_mols:
packing_molecule.merge(mol)
def merge_contacts(molecular_contacts):
for contact in molecular_contacts:
mol = contact.molecule2()
if mol not in core_mols:
packing_molecule.merge(mol)
for molecule_subset in core_mols:
if len(atom_selection) != 0:
# match the atoms to the molecule
atoms_to_test = []
for j in range(molecule_subset.natoms()):
# To replicate the comparison using numpy.allclose() we use,
# for all atom pairs:
# absolute(`a` - `b`) <= (`atol` + `rtol` * absolute(`b`))
# where rtol=1.e-5, atol=1.e-8
subset_atom = molecule_subset.atom(j)
o = subset_atom.site().orth()
for at in atom_selection:
if abs(at.coordinates.x - o.x()) <= 1.e-8 + 1.e-5 * abs(o.x()) and \
abs(at.coordinates.y - o.y()) <= 1.e-8 + 1.e-5 * abs(o.y()) and \
abs(at.coordinates.z - o.z()) <= 1.e-8 + 1.e-5 * abs(o.z()):
# atoms_to_test.append(Atom(_atom=molecule_subset.atom(j))._atom)
atoms_to_test.append(subset_atom)
break
if not atoms_to_test:
continue
# just get the contacts for the listed atoms
molecular_contacts = crystal_view.find_contacts(
molecule_subset, atoms_to_test, contact_criterion)
merge_contacts(molecular_contacts)
else:
# contacts for all components
molecular_contacts = crystal_view.find_contacts(
molecule_subset, contact_criterion)
merge_contacts(molecular_contacts)
return Molecule(self.identifier, packing_molecule).copy()
@property
def is_centrosymmetric(self):
'''Whether or not the crystal is centrosymmetric.'''
return self._crystal.cell().spacegroup().centrosymmetric()
@property
def is_sohncke(self):
'''Whether or not the crystal is in a Sohncke spacegroup.'''
return self._crystal.cell().spacegroup().sohncke()
[docs] def hbond_network(self, hbond_criterion=None, repetitions=1):
'''The molecule that results from expanding all hbonds in the crystal.
:param hbond_criterion: a :class:`ccdc.molecule.Molecule.HBondCriterion` instance.
:param repetitions: the number of times hbond contacts will be expanded.
'''
if hbond_criterion is None:
hbond_criterion = Molecule.HBondCriterion()
hbond_criterion.require_hydrogens = True
hbond_criterion.path_length_range = (3, 999)
csv = ChemistryLib.CrystalStructureView_instantiate(self._csv)
conv = ChemistryLib.CrystalStructureContactsView(csv)
conv.set_contact_criterion(hbond_criterion._hbc.hbond_criterion())
for i in range(repetitions):
contacts = conv.contacts()
conv.expand_contacts(contacts)
conv.delete_hanging_contacts()
mols = [csv.molecule(i) for i in range(csv.nmolecules())]
shell = Molecule(identifier=self.identifier)
for m in mols:
shell.add_molecule(Molecule(_molecule=m))
return shell
@staticmethod
def _decode_inclusion(inclusion):
if inclusion.lower().startswith('centroid'):
crit = ChemistryLib.AtomIfCentroidIncluded()
elif inclusion.lower().startswith('all'):
crit = ChemistryLib.AtomIfAllIncluded()
elif inclusion.lower().startswith('any'):
crit = ChemistryLib.AtomIfAnyIncluded()
elif inclusion.lower().startswith('only'):
crit = ChemistryLib.AtomOnlyIfIncluded()
else:
raise TypeError('Invalid inclusion criterion %s' % inclusion)
return crit
[docs] def packing(self, box_dimensions=((0, 0, 0), (1, 1, 1)), inclusion='CentroidIncluded'):
'''A molecule which fills some multiple of the unit cell of the crystal.
The atoms to include are specified with:
- 'CentroidIncluded' where whole molecules will be included if their centroid is within the box dimensions,
- 'AllAtomsIncluded' where whole molecules will be included only if all atoms of the molecule lie within the box,
- 'AnyAtomIncluded' where whole molecules will be included if any atom of the molecule lies within the box,
- 'OnlyAtomsIncluded' where all and only the atoms lying within the box will be included.
The default is to fill the unit cell.
:parameter box_dimensions: a pair of triples being the minimum and maximum multipliers of the unit cell axes to
fill with the crystal's molecules.
:parameter inclusion: one of 'CentroidIncluded', 'AllAtomsIncluded', 'AnyAtomIncluded', or 'OnlyAtomsIncluded'
where all and only the atoms lying within the box will be included.
:returns: a :class:`ccdc.molecule.Molecule` instance.
'''
csv = ChemistryLib.CrystalStructureView_instantiate(self._csv)
psv = ChemistryLib.CrystalStructurePackingView(csv)
if psv.packable():
csv.set_inclusion_condition(self._decode_inclusion(inclusion))
box = MathsLib.Box(box_dimensions[0], box_dimensions[1])
psv.set_packing_box(box)
psv.set_packing(True)
else:
raise RuntimeError('The crystal is not packable.')
mol = ChemistryLib.Molecule3d_instantiate()
for i in range(csv.nmolecules()):
mol.add(csv.molecule(i))
return Molecule(self.identifier, _molecule=mol)
[docs] def slicing(self, plane, thickness=10.0, width=20.0, displacement=0.0, inclusion='CentroidIncluded'):
'''The molecules of the crystal subtended by slicing along a plane.
:parameter plane: a slicing plane.
:parameter thickness: the thickness of the slicing view.
:parameter width: the width of the slicing view.
:parameter displacement: the displacement of the slicing view.
:parameter inclusion: one of 'CentroidIncluded', 'AllAtomsIncluded', 'AnyAtomIncluded', or 'OnlyAtomsIncluded'.
:returns: a :class:`ccdc.molecule.Molecule` instance.
'''
csv = ChemistryLib.CrystalStructureView_instantiate(self._csv)
csv.set_inclusion_condition(self._decode_inclusion(inclusion))
sv = ChemistryLib.CrystalStructureSlicingView(csv)
p = plane._plane
sv.set_plane(p)
sv.set_thickness(thickness)
sv.set_width(width)
sv.set_displacement(displacement)
sv.set_slicing(True)
mol = ChemistryLib.Molecule3d_instantiate()
for i in range(csv.nmolecules()):
mol.add(csv.molecule(i))
return Molecule(self.identifier, _molecule=mol)
[docs] def polymer_expansion(self, atoms=None, repetitions=1):
'''The molecule resulting from the expansion of polymeric bonds in the crystal.
:parameter atoms: an iterable of :class:`ccdc.molecule.Atom` instances whose polymeric bonds are to be expanded.
:parameter repetitions: the number of expansions to be performed.
:returns: :class:`ccdc.molecule.Molecule` instance.
'''
pe = ChemicalAnalysisLib.PolymerExpander()
pe.set_monomer(self._crystal)
if atoms is None:
for i in range(repetitions):
pe.expand_all()
else:
_m = self._crystal.editable_molecule()
for i in range(repetitions):
for a in atoms:
# find expandable bond containing a matching atom
polys = [_m.bond(i) for i in range(_m.nbonds()) if _m.bond(i).polymeric()]
for p in polys:
if p.atom1().label() == a.label:
# Expand this
pe.expand_from(p.atom1())
break
elif p.atom2().label() == a.label:
pe.expand_from(p.atom2())
break
mol = Molecule(self.identifier, _molecule=self._crystal.editable_molecule().clone())
pe.reset()
return mol
[docs] def to_string(self, format='mol2'):
'''Convert the crystal to a mol2 representation.
:param format: 'mol2', 'sdf', 'cif', 'mmcif' or 'smiles'
:returns: string representation in the appropriate format
:raises: TypeError if format is not as above.
'''
stream = FileFormatsLib.ostringstream()
mf = _file_object_factory(format)
mf.set(self._crystal, UtilitiesLib.DatabaseEntryIdentifier(self.identifier))
mf.write(stream)
return stream.str()
[docs] @staticmethod
def from_string(s, format=''):
'''Create a crystal from a string representation.
The format will be auto-detected if not specified.
:param s: string representation of a crystal or molecule
:param format: one of 'mol2', 'sdf', 'mol', 'cif', 'mmcif' or 'smiles'
:returns: a :class:`ccdc.crystal.Crystal`
:raises: TypeError if the format string is not '', 'mol2', 'sdf', 'mol', 'cif' or 'smiles'.
:raises: RuntimeError if the string representation is incorrectly formatted
'''
if format == '':
format = _detect_format(s)
if format in ['sdf', 'mol']:
if 'V2000' not in s and 'M END' not in s and '$$$$' not in s:
raise RuntimeError('This does not appear to be an SDF format string: %s' % s)
if format == 'cif':
format = 'cifreader'
stream = FileFormatsLib.istringstream(str(s))
mf = _file_object_factory(format)
mf.read(stream)
if format in ['mol', 'sdf']:
ident = mf.molecule_name()
else:
ident = mf.identifier().str()
try:
c = Crystal(mf.crystal_structure(), ident)
except RuntimeError as exc:
if 'There is no data' in str(exc):
c = Crystal(ChemistryLib.ConcreteCrystalStructure(), ident)
else:
raise RuntimeError('Invalid content for {} formatted string: {}'.format(format, s))
return c
@property
def inverted_structure(self):
'''Return an inverted crystal structure
:returns: a :class:`ccdc.crystal.Crystal` instance with inverted structure
:raises: RuntimeError if unable to invert crystal structure
'''
new_cs = self._crystal.clone()
if new_cs.invert_crystal_structure():
return type(self)(new_cs, self._identifier)
else:
raise RuntimeError('Unable to invert crystal structure')
[docs] def miller_indices(self, h, k, l):
'''
The :class:`ccdc.crystal.Crystal.MillerIndices` instance corresponding to integers h, k, l.
:param h: Miller index
:param k: Miller index
:param l: Miller index
:return: :class:`ccdc.crystal.Crystal.MillerIndices` instance.
:raises: TypeError if Miller indices are 0, 0, 0.
'''
return Crystal.MillerIndices(h, k, l, self)
[docs] def generate_inchi(self, include_stereo=True, add_hydrogens=True):
'''Return a :class:`ccdc.descriptors.MolecularDescriptors.InChIGenerator.InChI` object for the molecule
:param include_stereo: include stereochemistry (True, default) or ignore stereochemistry (False) when generating the InChI object
:param add_hydrogens: add hydrogens (True, default) or not add hydrogens (False) when generating the InChI object
:returns: a :class:`ccdc.descriptors.MolecularDescriptors.InChIGenerator.InChI` object
'''
from ccdc.descriptors import MolecularDescriptors
gen = MolecularDescriptors.InChIGenerator()
return gen.generate(self, include_stereo, add_hydrogens)
@property
def disorder(self):
'''The disorder object of the crystal.
Returns None if crystal structure has no disorder.
:returns: A :class:`ccdc.crystal.Crystal.Disorder` object.
'''
if self.has_disorder:
return Crystal.Disorder(self._csv)
else:
return None
##################################################################################
[docs]class PackingSimilarity(object):
'''Compare crystal packing similarities.
The crystal packing similarity feature is available only to CSD-Materials and CSD-Enterprise
users.
'''
_telemetry = 0
[docs] class Settings(object):
'''Settings for Packing similarity.'''
def __init__(self, _settings=None):
'''Initialise settings.'''
if _settings is None:
_settings = PackingSimilarityLib.PackingSimilaritySearchOptions()
self._settings = _settings
self.allow_molecular_differences = False
self.allow_artificial_inversion = True
@property
def angle_tolerance(self):
'''Maximum difference between reference and target angles.'''
return self._settings.angle_tolerance().degrees()
@angle_tolerance.setter
def angle_tolerance(self, value):
'''Set the angle tolerance.'''
self._settings.set_angle_tolerance(
MathsLib.Angle(value, MathsLib.Angle.DEGREES)
)
@property
def distance_tolerance(self):
'''Maximum difference between reference and target distances as a decimal fraction.'''
return self._settings.distance_tolerance()
@distance_tolerance.setter
def distance_tolerance(self, distance):
self._settings.set_distance_tolerance(distance)
@property
def allow_molecular_differences(self):
'''Whether different compounds may be compared.'''
return not self._settings.require_same_compound()
@allow_molecular_differences.setter
def allow_molecular_differences(self, value):
self._settings.set_require_same_compound(not value)
@property
def packing_shell_size(self):
'''The number of molecules in the packing shell.'''
return self._settings.packing_shell_size()
@packing_shell_size.setter
def packing_shell_size(self, value):
self._settings.set_packing_shell_size(value)
@property
def match_entire_packing_shell(self):
'''Whether or not all molecules of the shell have to be matched.'''
return self._settings.match_entire_packing_shell()
@match_entire_packing_shell.setter
def match_entire_packing_shell(self, value):
self._settings.set_match_entire_packing_shell(value)
@property
def ignore_hydrogen_positions(self):
'''Whether or not H positions should be ignored.'''
return self._settings.ignore_hydrogen_positions()
@ignore_hydrogen_positions.setter
def ignore_hydrogen_positions(self, value):
self._settings.set_ignore_hydrogen_positions(value)
@property
def ignore_bond_types(self):
'''Whether or not to take account of bond types.'''
return self._settings.ignore_bond_types()
@ignore_bond_types.setter
def ignore_bond_types(self, value):
self._settings.set_ignore_bond_types(value)
@property
def ignore_hydrogen_counts(self):
'''Whether or not to ignore hydrogen counts.'''
return self._settings.ignore_hydrogen_counts()
@ignore_hydrogen_counts.setter
def ignore_hydrogen_counts(self, value):
self._settings.set_ignore_hydrogen_counts(value)
@property
def ignore_bond_counts(self):
'''Whether or not to take account of bond counts.'''
return self._settings.ignore_bond_counts()
@ignore_bond_counts.setter
def ignore_bond_counts(self, value):
self._settings.set_ignore_bond_counts(value)
@property
def ignore_smallest_components(self):
'''Whether or not to take account of solvents.'''
return self._settings.ignore_smallest_components()
@ignore_smallest_components.setter
def ignore_smallest_components(self, value):
self._settings.set_ignore_smallest_components(value)
@property
def show_highest_similarity_result(self):
'''Control number of results.
For structures with Z' > 1 there will be more than one result.
'''
return self._settings.show_highest_similarity_result()
@show_highest_similarity_result.setter
def show_highest_similarity_result(self, value):
self._settings.set_show_highest_similarity_result(value)
@property
def skip_when_identifiers_equal(self):
'''Do not compare structures with the same identifier.'''
return self._settings.skip_when_identifiers_equal()
@skip_when_identifiers_equal.setter
def skip_when_identifiers_equal(self, value):
self._settings.set_skip_when_identifiers_equal(value)
@property
def molecular_similarity_threshold(self):
'''Do not compare structures whose similarity is lower than this value.'''
return self._settings.molecular_similarity_threshold()
@molecular_similarity_threshold.setter
def molecular_similarity_threshold(self, value):
self._settings.set_molecular_similarity_threshold(value)
@property
def allow_artificial_inversion(self):
'''Whether or not to invert a structure where there is no inversion symmetry.'''
return self._settings.force_artificial_inversion()
@allow_artificial_inversion.setter
def allow_artificial_inversion(self, value):
self._settings.set_force_artificial_inversion(value)
@property
def timeout_ms(self):
'''A timeout for the packing similarity calculation. 0 means no timeout.'''
return self._settings.timeout_ms()
@timeout_ms.setter
def timeout_ms(self, value):
self._settings.set_timeout_ms(value)
[docs] class Comparison(object):
'''The result of a comparison.'''
def __init__(self, _comp, reference, target):
'''Initialise the comparison.'''
self._comp = _comp
feature = self._comp.packing_shell_motif().crystal_feature()
feature_id = feature.source_crystal_identifier().str()
if feature_id == reference.identifier:
self.reference = reference
self.target = target
else:
self.reference = target
self.target = reference
@property
def rmsd(self):
'''RMSD of the reference and the target.'''
return self._comp.rmsd()
@property
def nmatched_molecules(self):
'''The number of molecules of the crystal matched.'''
return self._comp.nmatched_molecules()
@property
def packing_shell_size(self):
'''Number of molecules in the packing shell.'''
return self._comp.packing_shell_size()
[docs] def overlay_molecules(self):
'''Return the overlayed molecules.
:returns: a pair of packing shells, the first of which has been overlayed on the second
'''
best_match = self._comp.best_similarity_match()
feature = self._comp.packing_shell_motif().crystal_feature()
source_view = ChemistryLib.CrystalStructureView.instantiate(feature.source_crystal())
source = Molecule(feature.source_crystal_identifier().str())
mols = feature.feature()
for i, m in enumerate(mols):
if best_match.substructure_match(i).is_found():
at = m.atom(0)
whole = PackingSimilarityLib.find_molecule(at, source_view)
source.add_molecule(Molecule('', whole))
target_view = ChemistryLib.CrystalStructureView.instantiate(self.target._crystal)
mss = MotifSearchLib.MotifSearchStructure(target_view)
target = Molecule(self.target.identifier)
for i in range(best_match.nsubstructure_matches()):
match = best_match.substructure_match(i)
if match.is_found():
target.add_molecule(Molecule('', _molecule=mss.molecule(match)))
xform = self._comp.transformation()
mol_xform = ChemistryLib.MoleculeTransformation(target._molecule, xform)
mol_xform.apply()
return (source, target)
def __init__(self, settings=None):
'''Initialise the packing similarity instance.'''
if settings is None:
settings = PackingSimilarity.Settings()
self.settings = settings
if type(self)._telemetry == 0:
UtilitiesLib.ccdc_packing_similarity_telemetry()
type(self)._telemetry = 1
[docs] def compare(self, reference, target):
'''Compare two crystals.
:param reference: :class:`ccdc.crystal.Crystal`
:param target: :class:`ccdc.crystal.Crystal`
:returns: :class:`ccdc.crystal.PackingSimilarity.Comparison`,
a tuple of comparisons if there are more than one or ``None`` if no comparison was possible
'''
csv0 = ChemistryLib.CrystalStructureView.instantiate(reference._crystal)
csv1 = ChemistryLib.CrystalStructureView.instantiate(target._crystal)
msi0 = PackingSimilarityLib.MotifSurveyorCrystalItem(reference.identifier, csv0)
msi1 = PackingSimilarityLib.MotifSurveyorCrystalItem(target.identifier, csv1)
self._packing_similarity_search = PackingSimilarityLib.PackingSimilaritySearch()
self._packing_similarity_search.set_options(self.settings._settings)
self._packing_similarity_search.add_reference_structure(msi0)
self._packing_similarity_search.add_comparison_structure(msi1)
self._packing_similarity_search.construct_pairwise_comparisons()
l = self._packing_similarity_search.run_comparison(msi0, msi1)
if len(l) == 1:
return PackingSimilarity.Comparison(l[0], reference, target)
elif len(l) == 0:
return None
else:
return tuple(PackingSimilarity.Comparison(x, reference, target) for x in l)
@_if_pdd_is_licensed
def compare_with_filtering(self, reference_crystals, comparison_crystals,
atom_shell_size, pdd_threshold, amd_threshold,
apply_cell_expansion=False, cell_volume_expansion_threshold=12.5):
'''Compare two lists of crystals with filtering using PackingSimilarityGeometricFilter
:param reference_crystals: [:class:`ccdc.crystal.Crystal`]
:param comparison_crystals: [:class:`ccdc.crystal.Crystal`]
:param atom_shell_size: the integer specifying the number of nearest neigbhours
:param pdd_threshold: the positive float to indicate the threshold for PDD comparision
:param amd_threshold: the positive float to indicate the threshold for AMD comparision
:param apply_cell_expansion: the flag to indicate if to expand cell in the PDD calculation
:param cell_volume_expansion_threshold: the positive float to indicate the threshold for the cell expansion
:returns: [tuples of (ref_crystal, comp_crystal, :class:`ccdc.crystal.PackingSimilarity.Comparison`,
emd_distance, amd_distance, filtered flag)],
a list of tuples of (ref_crystal, comp_crystal, comparisons, emd_distance, amd_distance, filtered),
this list comprises all the pairs of ref_crystals and comp_crystals, and their EMD, AMD distances,
the filtered flag indicate if the pair has passed the filtering or not. The filtered out pairs
didn't go through compare(), thus the 'comparisons' in the list for such pairs are empty tuples
'''
settings = PackingSimilarityGeometricFilter.Settings()
settings.ignore_smallest_components = self.settings.ignore_smallest_components
settings.ignore_hydrogen_positions = self.settings.ignore_hydrogen_positions
settings.skip_when_identifiers_equal = self.settings.skip_when_identifiers_equal
settings.atom_shell_size = atom_shell_size
settings.pdd_threshold = pdd_threshold
settings.amd_threshold = amd_threshold
settings.apply_cell_expansion = apply_cell_expansion
settings.cell_volume_expansion_threshold = cell_volume_expansion_threshold
# filtering first
psf = PackingSimilarityGeometricFilter(settings)
comparison_maps = psf.filter_structures(reference_crystals, comparison_crystals)
# get the final list of results of (ref_crystal, comp_crystal, Comparison, emd, amd, filtered),
# some go through self.compare(), some just use empty tuple
compare_result = []
filtered_result = []
for comparison_map in comparison_maps:
for i in range(len(comparison_map.comp_structures)):
comp_structure = comparison_map.comp_structures[i]
if comparison_map.filtered[i] == 0:
compare_result.append((comparison_map.ref_structure, comp_structure,
self.compare(comparison_map.ref_structure.crystal, comp_structure.crystal),
comparison_map.emd_distances[i], comparison_map.amd_distances[i], False))
else:
filtered_result.append((comparison_map.ref_structure, comp_structure, tuple(),
comparison_map.emd_distances[i], comparison_map.amd_distances[i], True))
return compare_result + filtered_result
##################################################################################
@_if_pdd_is_licensed
class PackingSimilarityGeometricFilter(object):
'''Compare crystal packing similarities. This class does the filtering before going to PackingSimilarity
The crystal packing similarity feature is available only to CSD-Materials and CSD-Enterprise
users.
'''
_telemetry = 0
class Settings(object):
'''Settings for PDD(Pointwise Distance Distributions) and AMD(Average Minimum Distance).'''
def __init__(self, ignore_hydrogen_positions=True, skip_when_identifiers_equal=True,
ignore_smallest_components=True, atom_shell_size=100, pdd_threshold=0.35, amd_threshold=0.35,
apply_cell_expansion=False, cell_volume_expansion_threshold=12.5):
'''Initialise settings.'''
self._ignore_hydrogen_positions = ignore_hydrogen_positions
self._skip_when_identifiers_equal = skip_when_identifiers_equal
self.ignore_smallest_components = ignore_smallest_components
self._atom_shell_size = atom_shell_size
self._pdd_threshold = pdd_threshold
self._amd_threshold = amd_threshold
self._apply_cell_expansion = apply_cell_expansion
self._cell_volume_expansion_threshold = cell_volume_expansion_threshold
@property
def ignore_hydrogen_positions(self):
'''whether to ignore hydrogen positions.'''
return self._ignore_hydrogen_positions
@ignore_hydrogen_positions.setter
def ignore_hydrogen_positions(self, value):
'''Set whether to ignore hydrogen positions.'''
if type(value) is not bool:
raise TypeError("The value for ignore_hydrogen_positions need to be True or False!")
self._ignore_hydrogen_positions = value
@property
def ignore_smallest_components(self):
'''whether to ignore smallest components.'''
return self._ignore_smallest_components
@ignore_smallest_components.setter
def ignore_smallest_components(self, value):
'''Set whether to ignore smallest components.'''
if type(value) is not bool:
raise TypeError("The value for ignore_hydrogen_positions need to be True or False!")
self._ignore_smallest_components = value
@property
def skip_when_identifiers_equal(self):
'''whether to skip when identifiers equal.'''
return self._skip_when_identifiers_equal
@skip_when_identifiers_equal.setter
def skip_when_identifiers_equal(self, value):
'''Set whether to skip when identifiers equal.'''
if type(value) is not bool:
raise TypeError("The value for ignore_hydrogen_positions need to be True or False!")
self._skip_when_identifiers_equal = value
@property
def atom_shell_size(self):
'''the atom shell size for the PDD calculation'''
return self._atom_shell_size
@atom_shell_size.setter
def atom_shell_size(self, value):
'''Set the atom shell size for the PDD calculation'''
if type(value) is not int:
raise TypeError("The value for atom_shell_size need to be int between(0, 1000)!")
if value <= 0 or value >= 1000:
raise TypeError("The value for atom_shell_size need to be int between(0, 1000)!")
self._atom_shell_size = value
@property
def pdd_threshold(self):
'''the PDD threshold for the filtering comparison'''
return self._pdd_threshold
@pdd_threshold.setter
def pdd_threshold(self, value):
'''Set the PDD threshold for the filtering comparison'''
if type(value) is not float:
raise TypeError("The value for pdd_threshold need to be float!")
if value < 0.0:
raise TypeError("The value for pdd_threshold need to be positive!")
self._pdd_threshold = value
@property
def amd_threshold(self):
'''the AMD threshold for the filtering comparison'''
return self._amd_threshold
@amd_threshold.setter
def amd_threshold(self, value):
'''Set the AMD threshold for the filtering comparison'''
if type(value) is not float:
raise TypeError("The value for amd_threshold need to be float!")
if value < 0.0:
raise TypeError("The value for amd_threshold need to be positive!")
self._amd_threshold = value
@property
def apply_cell_expansion(self):
'''whether to apply cell expansion.'''
return self._apply_cell_expansion
@apply_cell_expansion.setter
def apply_cell_expansion(self, value):
'''Set whether to apply cell expansion.'''
if type(value) is not bool:
raise TypeError("The value for apply_cell_expansion need to be True or False!")
self._apply_cell_expansion = value
@property
def cell_volume_expansion_threshold(self):
'''The threshold for cell volume expansion'''
return self._cell_volume_expansion_threshold
@cell_volume_expansion_threshold.setter
def cell_volume_expansion_threshold(self, value):
'''Set the threshold for cell volume expansion'''
value = float(value)
if type(value) is not float:
raise TypeError("The value for cell_volume_expansion_threshold need to be a float between(0.0, 100.0)!")
if value <= 0.0 or value >= 100.0:
raise TypeError("The value for cell_volume_expansion_threshold need to be a float between(0.0, 100.0)!")
self._cell_volume_expansion_threshold = value
_settings = Settings()
class PeriodicSet(object):
def __init__(self, settings=None):
'''Initialise PeriodicSet with settings
use PackingSimilarityGeometricFilter._settings if settins is None'''
if settings is None:
settings = PackingSimilarityGeometricFilter._settings
self._settings = settings
self._ps = None
def from_crystal(self, crystal, name):
''' Get the PeriodicSet from crystal
>>> from ccdc.io import EntryReader
>>> csd_reader = EntryReader('CSD')
>>> ibprac01 = csd_reader.crystal('IBPRAC01')
>>> ps = PackingSimilarityGeometricFilter().PeriodicSet()
>>> ps.from_crystal(ibprac01, 'IBPRAC01')
>>> print(ps._ps.identifier())
IBPRAC01
'''
self._ps = PddLib.PeriodicSetCalculator.from_ccdc_crystal(
crystal._crystal, name, self._settings.ignore_hydrogen_positions, self._settings.ignore_smallest_components)
def from_entry(self, entry):
''' Get the PeriodicSet from entry
>>> from ccdc.io import EntryReader
>>> csd_reader = EntryReader('CSD')
>>> qaxmeh = csd_reader.entry('QAXMEH')
>>> ps = PackingSimilarityGeometricFilter().PeriodicSet()
>>> ps.from_entry(qaxmeh)
>>> print(ps._ps.identifier())
QAXMEH
'''
self._ps = PddLib.PeriodicSetCalculator.from_ccdc_entry(
entry._entry, self._settings.ignore_hydrogen_positions, self._settings.ignore_smallest_components)
class PDD(object):
''' PointwiseDistanceDistribution '''
def __init__(self, periodic_set=None, settings=None):
'''Initialise PDD with PeriodicSet and settings,
use PackingSimilarityGeometricFilter._settings if settins is None'''
if settings is None:
settings = PackingSimilarityGeometricFilter._settings
self._settings = settings
if periodic_set is None:
periodic_set = PackingSimilarityGeometricFilter.PeriodicSet()
self._periodic_set = periodic_set
self._pdd = MathsLib.PointwiseDistanceDistribution(self._periodic_set._ps, self._settings.atom_shell_size)
@property
def matrix(self):
''' Return the PDD matrix property. '''
return self._pdd.matrix()
def distance(self, pdd2):
''' Return the EMD distance. '''
return self._pdd.distance(pdd2._pdd)
class AMD(object):
'''AverageMinimumDistance'''
def __init__(self, periodic_set=None, settings=None):
'''Initialise the AMD with PeriodicSet and settings,
use PackingSimilarityGeometricFilter._settings if settins is None'''
if settings is None:
settings = PackingSimilarityGeometricFilter._settings
self._settings = settings
if periodic_set is None:
periodic_set = PackingSimilarityGeometricFilter.PeriodicSet()
self._periodic_set = periodic_set
self._amd = MathsLib.AverageMinimumDistance(self._periodic_set._ps, self._settings.atom_shell_size)
def from_pdd(self, pdd):
'''Reset the AMD from PDD'''
self._amd = MathsLib.AverageMinimumDistance(pdd.matrix)
def from_ps(self, ps):
'''Reset the AMD from a PeriodicSet'''
self._amd = MathsLib.AverageMinimumDistance(ps._ps, self._settings.atom_shell_size)
@property
def vector(self):
'''Return the AMD vector property'''
return self._amd.to_vector()
def distance(self, amd=None):
'''Return the distance of this AMD to another'''
if amd is None:
amd = PackingSimilarityGeometricFilter.AMD()
return self._amd.distance(amd._amd)
class CrystalStructureIdentifier:
'''The structure that combines the crystal structure and its identifier'''
def __init__(self, cs_id):
self._cs = Crystal(cs_id.crystal_structure(), cs_id.identifier())
@property
def crystal(self):
return self._cs
@property
def identifier(self):
return self._cs._identifier
class ComparisonMap:
'''The structure that maps the reference structure to the comparison structures
with the amd and emd distances and the flag if the pair has passed the filtering '''
def __init__(self, filtered_structure):
self._ref_structure = filtered_structure.reference_structure()
self._comp_structures = filtered_structure.comparison_structures()
self._emd_distances = filtered_structure.emd()
self._amd_distances = filtered_structure.amd()
self._filtered = filtered_structure.filtered_out()
@property
def ref_structure(self):
return PackingSimilarityGeometricFilter.CrystalStructureIdentifier(self._ref_structure)
@property
def comp_structures(self):
return [PackingSimilarityGeometricFilter.CrystalStructureIdentifier(structure)
for structure in self._comp_structures]
@property
def emd_distances(self):
return [emd_dist for emd_dist in self._emd_distances]
@property
def amd_distances(self):
return [amd_dist for amd_dist in self._amd_distances]
@property
def filtered(self):
return [is_filtered for is_filtered in self._filtered]
def __init__(self, settings=None):
'''Initialise the PackingSimilarityGeometricFilter instance with settings.'''
if settings is None:
settings = PackingSimilarityGeometricFilter.Settings()
PackingSimilarityGeometricFilter._settings = settings
def filter_structures(self, reference_crystals, comparison_crystals):
ref_crystals = [ref_cs._crystal for ref_cs in reference_crystals]
ref_ids = [ref_cs.identifier for ref_cs in reference_crystals]
comp_crystals = [comp_cs._crystal for comp_cs in comparison_crystals]
comp_ids = [comp_cs.identifier for comp_cs in comparison_crystals]
maps = PddLib.PackingSimilarityGeometricFilter.filter_structures(
ref_crystals, ref_ids, comp_crystals, comp_ids, self._settings.amd_threshold,
self._settings.pdd_threshold, self._settings.ignore_hydrogen_positions,
self._settings.ignore_smallest_components, self._settings.skip_when_identifiers_equal,
self._settings.atom_shell_size, self._settings.apply_cell_expansion,
self._settings.cell_volume_expansion_threshold)
comparison_maps = [PackingSimilarityGeometricFilter.ComparisonMap(map) for map in maps]
return comparison_maps
##################################################################################
