Source code for ccdc.crystal

#
# 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 Contact(Molecule.Contact): '''A crystallographic contact.''' _type_map = dict(( (ChemistryLib.MolecularContact.UNKNOWN_CONTACT, "Unknown"), (ChemistryLib.MolecularContact.VDW_CONTACT, "Van der Waals"), (ChemistryLib.MolecularContact.DISTANCE_CONTACT, "Distance"), (ChemistryLib.MolecularContact.HBOND_CONTACT, "HBond"), (ChemistryLib.MolecularContact.REAL_BOND_CONTACT, "Real Bond"), (ChemistryLib.MolecularContact.CUSTOM_CONTACT, "Custom"), )) 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 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)
[docs] def contacts(self, intermolecular='Any', distance_range=(-5.0, 0.0), path_length_range=(4, 999)): '''The collection of crystallographic contacts in this crystal structure. A crystallographic contact may span two symmetry related molecules of the crystal. Distance limits for nonbonded contacts are defined as distances relative to the sum of van der Waals Radii, *i.e.*: v(X) + v(Y) - p < d(X...Y) < v(X) + v(Y) - q where *v(X)* and *v(Y)* are the vdW radii of *X* and *Y* respectively and *p* and *q* are user-specified values in Angstroms (*e.g.*, if *p* = 0.5 and *q* = 0.1, the contact must be at least 0.1 Angstroms shorter than the sum of vdW radii but no more than 0.5 Angstroms shorter). Standard van der Waals radii are assigned to the common elements. They are taken from Bondi, *J.Phys.Chem.*, **68**, 441, 1964. Other elements are assigned van der Waals radii of 2.0 Angstroms. :param intermolecular: 'Intramolecular', 'Intermolecular' or 'Any'. :param distance_range: Minimum and maximum distances considered acceptable for a contact to be formed (Angstroms). :param path_length_range: Minimum and maximum values for length of path between the contact atoms :returns: a tuple of :class:`ccdc.crystal.Crystal.Contact` ''' if intermolecular.lower().startswith('intra'): inter = ChemistryLib.ContactCriterion.INTRAMOLECULAR elif intermolecular.lower().startswith('inter'): inter = ChemistryLib.ContactCriterion.INTERMOLECULAR else: return self.contacts('inter', distance_range, path_length_range) + \ self.contacts('intra', distance_range, path_length_range) hbc = ChemistryLib.VdWCriterion( max(distance_range), min(distance_range), 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.Contact(_view=view, _contact=contact) for h in ret: if ( _base_atom(hb.atoms[0]) == _base_atom(h.atoms[0]) and _base_atom(hb.atoms[-1]) == _base_atom(h.atoms[-1]) and round(hb.length, 5) == round(h.length, 5) # and #hb.symmetry_operators[0] == h.symmetry_operators[-1] and #hb.symmetry_operators[-1] == h.symmetry_operators[0] ): break if ( _base_atom(hb.atoms[0]) == _base_atom(h.atoms[-1]) and _base_atom(hb.atoms[-1]) == _base_atom(h.atoms[0]) and round(hb.length, 5) == round(h.length, 5) ): 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 contact_network(self, intermolecular='Any', distance_range=(-5.0, 0.0), vdw_corrected=True, path_length_range=(4, 999), repetitions=1): '''The molecule resulting from expanding all crystallographic contacts of the crystal. A crystallographic contact may span two symmetry related molecules of the crystal. :param intermolecular: 'Intramolecular', 'Intermolecular' or 'Any'. :param distance_range: Minimum and maximum distances considered acceptable for a contact to be formed. :param vdw_corrected: Whether the distances are relative to the Van der Waals radius of the atoms. :param path_length_range: Minimum and maximum values for length of path between the contact atoms :returns: a tuple of :class:`ccdc.crystal.Crystal.Contact` ''' if intermolecular.lower().startswith('intra'): inter = ChemistryLib.ContactCriterion.INTRAMOLECULAR elif intermolecular.lower().startswith('inter'): inter = ChemistryLib.ContactCriterion.INTERMOLECULAR else: inter = ChemistryLib.ContactCriterion.ANY hbc = ChemistryLib.VdWCriterion( max(distance_range), min(distance_range), inter, min(path_length_range), max(path_length_range) ) csv = ChemistryLib.CrystalStructureView_instantiate(self._crystal) conv = ChemistryLib.CrystalStructureContactsView(csv) conv.set_contact_criterion(hbc) 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: to_add = Molecule(_molecule=m) shell.add_molecule(to_add) return shell
[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)
##################################################################################