#
# 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')
##################################################################################
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)
def __eq__(self, other):
'''Equality of crystals.'''
if not isinstance(other, Crystal):
return False
if self.identifier != other.identifier:
return False
if self._crystal == other._crystal:
return True
result = self._crystal.compare(other._crystal, 0)
return result == 0
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)
csv.all_disorder_grouping()
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.3
: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.CrystalBondTypesAssigner()
assigner.assign(self._crystal)
[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()
elif inclusion.lower().startswith('unique'):
crit = ChemistryLib.AtomIfHasAnyUniqueSupercellPosition()
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,
- 'UniqueIncluded' where whole molecules will be included if their centroid is within the box dimensions and they contribute unique box cell positions.
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', 'OnlyAtomsIncluded', or 'UniqueIncluded'
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', 'OnlyAtomsIncluded', or 'UniqueIncluded'.
: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
[docs]
def copy(self):
'''Return a deep copy of the crystal.'''
return type(self)(self._crystal.clone(), self._identifier)
@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]
@staticmethod
def generate_reduced_crystal(crystal):
'''This generates a copy of crystal in Niggli reduced cell form.
:param crystal: input :class:`ccdc.crystal.Crystal`
:returns: :class:`ccdc.crystal.Crystal` in Niggli reduced cell form
'''
reduced_crystal = crystal.copy()
ChemicalAnalysisLib.ReducedCellTransformer.transform_to_reduced_cell(reduced_crystal._crystal)
return reduced_crystal
[docs]
def reduce_symmetry(self):
'''This generates a copy of the crystal in reduced symmetry form, that is a
reduced spacegroup symmetry form so that no molecules occupy special positions
:param crystal: input :class:`ccdc.crystal.Crystal`
:returns: :class:`ccdc.crystal.Crystal` in reduced symmetry form
'''
reduced_crystal = self.copy()
ChemicalAnalysisLib.CrystalStructureSymmetryReducer.reduce_spacegroup_symmetry(reduced_crystal._crystal)
return reduced_crystal
[docs]
def reduce_symmetry_to_p1(self):
'''This generates a copy of the crystal in reduced symmetry form.
The symmetry is reduced to p1 in a single step.
:param crystal: input :class:`ccdc.crystal.Crystal`
:returns: :class:`ccdc.crystal.Crystal` in reduced symmetry form
'''
reduced_crystal = self.copy()
ChemicalAnalysisLib.CrystalStructureSymmetryReducer.reduce_spacegroup_to_p1(reduced_crystal._crystal)
return reduced_crystal
##################################################################################
[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)