#
# 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 operator
import sys
from ccdc.molecule import Molecule, _file_object_factory
from ccdc.utilities import nested_class, _detect_format
from ccdc.utilities import _private_importer
with _private_importer():
import UtilitiesLib
import MathsLib
import ChemistryLib
import ChemicalAnalysisLib
import FileFormatsLib
import MotifSearchLib
import PackingSimilarityLib
##################################################################################
CellLengths = collections.namedtuple('CellLengths', ['a', 'b', 'c'])
CellAngles = collections.namedtuple('CellAngles', ['alpha', 'beta', 'gamma'])
[docs]class Crystal(object):
'''Represents a crystal.'''
[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
def __init__(self, crystal, identifier):
'''Create crystal using a ChemistryLib.CrystalStructure and an identifier.'''
self._crystal = crystal
self._identifier = identifier
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 reduced_cell(self):
'''The reduced cell of the crystal.'''
return Crystal.ReducedCell(self._crystal.cell())
@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.
'''
m = self.disordered_molecule
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.
'''
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._crystal)
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 = ChemistryLib.CrystalStructureView_instantiate(self._crystal)
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
'''
vf = ChemicalAnalysisLib.CrystalVoidsFinder()
atoms = vf.calculate_atoms(self._crystal)
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._crystal, 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 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=[]):
"""
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.
:returns: :class:`ccdc.molecule.Molecule` containing the atoms within
the distance cutoff excluding any atoms from the molecules
defined in the atom selection.
"""
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()
for i in range( crystal_view.nmolecules() ):
molecule_subset = crystal_view.molecule(i)
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)
for contact in molecular_contacts:
packing_molecule.merge( contact.molecule2() )
else:
# contacts for all components
molecular_contacts = crystal_view.find_contacts(
molecule_subset, contact_criterion)
for contact in molecular_contacts:
packing_molecule.merge( contact.molecule2() )
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._crystal)
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._crystal)
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._crystal)
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' or 'cif'
: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' or 'cif'
: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 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]class PackingSimilarity(object):
'''Compare crystal packing similarities.
The crystal packing similarity feature is available only to CSD-Materials and CSD-Enterprise
users.
'''
[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)
[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
[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)
##################################################################################
