Source code for ccdc.particle

#
# This code is Copyright (C) 2020 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 :mod:`ccdc.particle` contains classes for particle surface analysis and slip planes calculation.

The :mod:`ccdc.particle` module contains two classes
:class:`ccdc.particle.Surface` and :class:`ccdc.particle.SlipPlanes` .

The :class:`ccdc.particle.Surface` class contains attributes relating
to the chemistry and topology of a particle surface.

.. code-block:: python

    >>> from ccdc import io
    >>> from ccdc.particle import Surface
    >>> crystal = io.CrystalReader('CSD').crystal('DEBXIT')
    >>> particle_surface = Surface(crystal, (1,0,0))
    >>> len(particle_surface.slab.atoms)
    414
    >>> particle_surface.slab.atoms[0].is_acceptor
    False
    >>> particle_surface.miller_indices.plane
    Plane(normal=Vector(1.000, -0.000, -0.000), distance=22.512)
    >>> round(len(list(particle_surface.topology.nodes)), -3)
    4000
    >>> tuple(round(x,0) for x in list(particle_surface.topology.nodes)[0])
    (22.0, 0.0, -0.0)
    >>> round(len(list(particle_surface.topology.triangles)), -3)
    7000
    >>> tuple(round(x,0) for x in list(particle_surface.topology.triangles)[0][0])
    (22.0, 0.0, -0.0)
    >>> particle_surface.descriptors.projected_area
    165.167
    >>> particle_surface.descriptors.rugosity
    1.372
    >>> particle_surface.descriptors.rmsd
    1.656
    >>> particle_surface.descriptors.skewness
    0.808
    >>> particle_surface.descriptors.kurtosis
    2.816
    >>> particle_surface.descriptors.hb_acceptors
    0.024
    >>> particle_surface.descriptors.hb_donors
    0.048
    >>> particle_surface.descriptors.unsatisfied_hb_donors
    0.012
    >>> particle_surface.descriptors.aromatic_bonds
    0.145

The :class:`ccdc.particle.SlipPlanes` class contains attributes relating
to the position of possible slip planes within a unit cell.

.. code-block:: python

    >>> from ccdc import io
    >>> from ccdc.particle import SlipPlanes
    >>> crystal = io.EntryReader('CSD').crystal("SCCHRN")
    >>> planes = SlipPlanes(crystal)
    >>> len(planes)
    1
    >>> round(planes[0].slab_separation, 3)
    1.095
    >>> [round(plane.slab_separation, 3) for plane in planes]
    [1.095]
    >>> planes.set_separation_threshold(-2)
    >>> [round(plane.slab_separation, 3) for plane in planes]
    [1.095, -0.063, -0.219, -0.34, -0.348, -0.348, -0.677, -1.232, -1.528, -1.724, -1.801, -1.867, -1.928]
    >>> planes[0].orientation
    MillerIndices(1, 0, 0)
    >>> planes[0].h_bonded
    False
    >>> planes[0].mid_plane
    Plane(normal=Vector(0.971, -0.000, 0.239), distance=9.285)
    >>> round(planes[0].offset, 3) + 0
    0.0
    >>> round(planes[0].repeat_distance, 3)
    9.285
'''
###############################################################################

import warnings

from ccdc import molecule
from ccdc.descriptors import GeometricDescriptors

from ccdc.utilities import nested_class, _private_importer

with _private_importer() as pi:
    pi.import_ccdc_module('ParticleLib')
    pi.import_ccdc_module('ChemistryLib')


###############################################################################


def no_coordinates(crystal):
    """
    :param crystal: Input crystal object
    :return True if crystal contains no atoms with coordinates, False otherwise
    """
    # Only atoms with coordinates are included in a Crystal object's molecule
    return len(crystal.molecule.atoms) == 0


def unknown_spacegroup(crystal):
    """
    :param crystal: Input crystal object
    :return True if spacegroup symbol is 'Unknown', False otherwise
    """
    return crystal.spacegroup_symbol == 'Unknown'


def check_for_siteless_atoms_or_unknown_bond_types(crystal, warnings, what):
    """
    :param crystal:  Input crystal object
    :param warnings: warnings object to add messages to
    :param what: string indicating whether it's for surface generation or slip plane generation
    """
    has_siteless_atoms = not crystal.molecule.all_atoms_have_sites
    has_unknown_bond_types = crystal.molecule.contains_unknown_bonds

    warning_text = ("{} - the {} "
                    "generated for this structure may not be reliable. We recommend "
                    "that you edit the structure to define atom sites and unknown "
                    "bond types prior to generation.").format(
        ParticleLib.basic_warning_text(has_siteless_atoms, has_unknown_bond_types), what)

    if has_siteless_atoms or has_unknown_bond_types:
        warnings.warn(warning_text, RuntimeWarning)


[docs]class Surface(): '''A class to represent a particle surface. Create a surface from a crystal and the given Miller indices or SlipPlane. :param crystal: The crystal whose surface we want to generate :type crystal: class:`ccdc.crystal.Crystal` :param miller_indices: Miller indices as a tuple of three integers from -9 to 9 :type miller_indices: tuple :param thickness_factor: slab thickness multiplier, defaults to automatic :type thickness_factor: float, optional :param surface_size: multiplier for u and v vectors, defaults to (1,1) :type surface_size: tuple, optional :param probe_radius: the probe radius for surface calculation, defaults to 1.2 :type probe_radius: float, optional :param grid_spacing: the grid spacing for surface calculation, defaults to 0.3 :type grid_spacing: float, optional :param offset: attempt to create a surface displaced by this distance from the default surface (offset=0) :type offset: float, optional :param slip_plane: A pre-calculated slip plane defining where to generate the surface :type slip_plane: class:`ccdc.particle.SlipPlane` '''
[docs] @nested_class('Surface') class Topology(): '''The topology of the surface, given by a collection of nodes in triangular mesh. ''' @staticmethod def __node_to_xyz(node): return (node.x(), node.y(), node.z()) def __init__(self, topology): self._topology = topology self._nodes = [] for nn in range(self._topology.nnodes()): self._nodes.append(Surface.Topology.__node_to_xyz( self._topology.node(nn))) self._triangles = [] for nt in range(self._topology.ntriangles()): triangle = self._topology.triangle(nt) node1 = Surface.Topology.__node_to_xyz(triangle.node1()) node2 = Surface.Topology.__node_to_xyz(triangle.node2()) node3 = Surface.Topology.__node_to_xyz(triangle.node3()) self._triangles.append((node1, node2, node3)) @property def nodes(self): '''Generator of nodes in the surface. A node is a tuple of three floats. ''' return (nn for nn in self._nodes) @property def triangles(self): '''Generator of triangles in the surface. A triangle is a tuple of three nodes. ''' return (nt for nt in self._triangles)
[docs] @nested_class('Surface') class Descriptors(): '''Descriptors of a surface.''' def __init__(self, surface, csv, cssv, contacts): self._surface_area = ParticleLib.surface_area(surface) self._projected_area = ParticleLib.projected_area(surface, cssv) self._rugosity = self._surface_area / self._projected_area density_info = ParticleLib.get_density_info( cssv, contacts, self._projected_area) self._hb_acceptors = density_info.acceptors_density self._hb_donors = density_info.donors_density self._unsatisfied_hb_donors = density_info.unsatisfied_donors_density self._aromatic_bonds = density_info.aromatic_bonds_density self._rmsd = ParticleLib.get_rmsd(surface, cssv) self._skewness = ParticleLib.get_skewness(surface, cssv, self._rmsd) self._kurtosis = ParticleLib.get_kurtosis(surface, cssv, self._rmsd) @property def hb_acceptors(self): '''Hydrogen bond acceptor density (count per Angstrom squared).''' return round(self._hb_acceptors, 3) @property def hb_donors(self): '''Hydrogen bond donor density (count per Angstrom squared).''' return round(self._hb_donors, 3) @property def unsatisfied_hb_donors(self): '''Unsatisfied Hydrogen bond donor density (count per Angstrom squared). Where an unsatisfied hydrogen bond donor is defined as a donor with a hydrogen atom that is surface terminating and not forming a hydrogen bond to the slab.''' return round(self._unsatisfied_hb_donors, 3) @property def aromatic_bonds(self): '''Aromatic bond density (count per Angstrom squared).''' return round(self._aromatic_bonds, 3) @property def surface_area(self): '''Topological surface area (in Angstrom squared).''' return round(self._surface_area, 3) @property def projected_area(self): '''Reticular area of the surface (in Angstrom squared).''' return round(self._projected_area, 3) @property def rugosity(self): '''Ratio between surface area and projected area.''' return round(self._rugosity, 3) @property def rmsd(self): '''Root mean squared deviation of the surface (in Angstrom)''' return round(self._rmsd, 3) @property def skewness(self): '''Measure of height deviation symmetry, where 0 is a symmetric distribution of peaks and valleys.''' return round(self._skewness, 3) @property def kurtosis(self): '''Pearson Kurtosis - a measure of the flatness of distribution.''' return round(self._kurtosis, 3)
[docs] @nested_class('Surface') class AtomProperties(): '''Atom properties of a surface.''' def __init__(self, cssv, contacts): self._a, self._d, self._ud = ParticleLib.get_atom_property_list(cssv, contacts) @property def hb_acceptors(self): '''List of atoms that are hydrogen bond acceptors''' return [molecule.Atom(_atom=a) for a in self._a] @property def hb_donors(self): '''List of atoms that are hydrogen bond donors''' return [molecule.Atom(_atom=a) for a in self._d] @property def unsatisfied_hb_donors(self): '''List of atoms that are unsatisfied hydrogen bond donors''' return [molecule.Atom(_atom=a) for a in self._ud]
def __init__(self, crystal, miller_indices=None, thickness_factor=0, surface_size=(1, 1), probe_radius=1.2, grid_spacing=0.3, offset=0.0, slip_plane=None): if unknown_spacegroup(crystal): raise RuntimeError("Insufficient information for surface analysis") if slip_plane and miller_indices: raise RuntimeError("Cannot specify both miller indices and a slip plane in the Surface constructor.") if slip_plane and offset: raise RuntimeError("Cannot specify both an offset and a slip plane in the Surface constructor.") check_for_siteless_atoms_or_unknown_bond_types(crystal, warnings, "surfaces") '''Construct a surface.''' self.crystal = crystal self.identifier = self.crystal.identifier self.surface_size = surface_size self.grid_spacing = grid_spacing self._offset = offset if slip_plane: self._miller_indices = slip_plane.orientation offset = slip_plane.offset else: if any([idx < -9 or idx > 9 for idx in miller_indices]): raise ValueError( "Miller indices must be integers between -9 and 9.") self._miller_indices = self.crystal.miller_indices(*miller_indices) d = self.miller_indices.plane.distance if offset < -0.99*d or offset > d: raise RuntimeError(f'offset out of range {round(-0.99*d,3),round(d,3)}') if any([mul not in [1, 2, 3, 4] for mul in surface_size]): raise ValueError( "Vector multipliers must be integers between 1 and 4.") csv = ChemistryLib.CrystalStructureView_instantiate( self.crystal._crystal) self._cssv = ParticleLib.CrystalStructureSurfaceView(csv) if thickness_factor == 0: thickness_factor = self._cssv.default_slab_thickness_factor(self._miller_indices._miller_indices) self.thickness_factor = thickness_factor self._cssv.set_slab_thickness_factor(thickness_factor) self._cssv.set_vector_u_multiplier(surface_size[0]) self._cssv.set_vector_v_multiplier(surface_size[1]) # periodic_view: true, preview: false, polymeric: false self._cssv.set_slicing_plane(self._miller_indices._miller_indices, True, False, False, offset) _topology, above_surface = ParticleLib.surface( self._cssv, grid_spacing, probe_radius) self._topology = Surface.Topology(_topology) self._above_surface = above_surface self._contacts = ParticleLib.get_surface_contacts(_topology, self._cssv, grid_spacing) self._atom_contacts = None self._surface_node_contacts = None self._descriptors = Surface.Descriptors(_topology, csv, self._cssv, self._contacts) self.atom_properties = Surface.AtomProperties(self._cssv, self._contacts) @property def surface_atoms(self): '''Returns a list of atoms found within VDW plus grid_spacing of the topology.''' return [molecule.Atom(_atom=a) for a in self._contacts.surface_atoms_] @property def periodic_vectors(self): ''' Returns the UVW vectors for constructing the periodic surface.''' u, v, w = [], [], [] for i in range(3): u.append(self._cssv.periodic_bounding_box().edge1()[i] / self.surface_size[0]) v.append(self._cssv.periodic_bounding_box().edge2()[i] / self.surface_size[1]) w.append(self._cssv.periodic_bounding_box().edge3()[i] / self.thickness_factor) u = GeometricDescriptors.Vector(*u) v = GeometricDescriptors.Vector(*v) w = GeometricDescriptors.Vector(*w) return [u, v, w] @property def topology(self): '''A :class:`ccdc.particle.Surface.Topology` object.''' return self._topology @property def miller_indices(self): '''The Miller plane indices used to define the particle Surface.''' return self._miller_indices @property def offset(self): '''The offset is the distance translated along the surface normal (in Angstroms).''' return self._offset @property def slicing_plane_distance(self): '''Distance from origin of reference plane (hkl and offset, or, SlipPlane) defining the Surface (in Angstroms).''' return self._cssv.get_slicing_plane().distance @property def slab(self): '''The particle slab represented as a single :class:`ccdc.molecule.Molecule`.''' surface_csv = self._cssv.periodic_crystal_view() self._slab = ChemistryLib.Molecule3d_instantiate() for i in range(surface_csv.nmolecules()): self._slab.add(surface_csv.molecule(i)) return molecule.Molecule(self.identifier, _molecule=self._slab) @property def slab_molecules(self): '''The particle slab represented as a list of the constituent Molecules.''' surface_csv = self._cssv.periodic_crystal_view() slab_molecules = [] for i in range(surface_csv.nmolecules()): mol = ChemistryLib.Molecule3d_instantiate(surface_csv.molecule(i)) slab_molecules.append(molecule.Molecule(self.identifier + str(i + 1), _molecule=mol)) return slab_molecules @property def descriptors(self): '''A :class:`ccdc.particle.Surface.Descriptors` object.''' return self._descriptors @property def atom_surface_node_contacts(self): '''Return a dict mapping :class:`ccdc.molecule.Atom` to a list of contacting surface node positions.''' if self._atom_contacts is None: self._atom_contacts = dict() for a in self._contacts.surface_atoms_: nodes = tuple(Surface.Topology._Topology__node_to_xyz(node) for node in ParticleLib.atom_node_contacts(self._contacts, a)) self._atom_contacts[molecule.Atom(_atom=a)] = nodes return self._atom_contacts @property def surface_node_atom_contacts(self): '''Return a dict mapping surface node positions to a list of contacting :class:`ccdc.molecule.Atom`.''' if self._surface_node_contacts is None: self._surface_node_contacts = dict() for node in ParticleLib.contact_nodes(self._contacts): atoms = tuple(molecule.Atom(_atom=a) for a in ParticleLib.node_atom_contacts(self._contacts, node)) self._surface_node_contacts[Surface.Topology._Topology__node_to_xyz(node)] = atoms return self._surface_node_contacts @property def _interaction_volume(self): '''Get the interaction volume for the surface (in Angstrom cubed).''' return ParticleLib.surface_interaction_volume(self._topology._topology, self._cssv, self._above_surface)
[docs]class SlipPlanes(): '''A list-like class to represent the slip planes in a crystal. The individual slip planes can be accessed by list-like index, slice or iteration. The initial separation threshold is 0.0. :param crystal: The crystal whose slip planes we want to calculate. :type crystal: class:`ccdc.crystal.Crystal` '''
[docs] class SlipPlane(): '''A class to represent a slip plane in a crystal. ''' def __init__(self, parent_planes, plane): self._planes = parent_planes self._plane = plane def __eq__(self, other): return str(self.mid_plane) == str(other.mid_plane) @property def slab_separation(self): '''Returns the shortest distance (in Angstroms) between the slabs of molecules on either side of the slip plane. A negative value indicates an overlap of the two slabs of molecules.''' return self._plane.slab_separation() @property def offset(self): '''Returns the perpendicular distance from this slip plane to the parallel Miller plane (in Angstroms).''' return self._plane.miller_plane_to_mid_plane_offset() @property def h_bonded(self): '''Returns whether or not intermolecular hydrogen bonds cross the slip plane.''' return self._plane.hbonds_between_layers() @property def orientation(self): '''Returns the indices of the Miller plane that sits parallel to the slip plane.''' indices = self._plane.miller_indices() return self._planes._crystal.miller_indices(indices.h(), indices.k(), indices.l()) @property def mid_plane(self): '''Returns the slip plane as a :class:`ccdc.descriptors.GeometricDescriptors.Plane`.''' mid = self._plane.mid_plane() return GeometricDescriptors.Plane(vector=None, distance=None, _plane=mid) @property def perpendicular_slip_planes(self): '''Returns a list of slip planes with orientations perpendicular to this one.''' perps = self._planes._planes_results.perpendicular_slip_planes(self._plane) return [SlipPlanes.SlipPlane(self._planes, perp) for perp in perps] @property def repeat_distance(self): '''Returns the perpendicular distance to the next equivalent slip plane (in Angstroms).''' return self._plane.repeat_distance() # define a property factory class member so that the properties can be declared once only # and the property names are available to SlipPlanes even if no SlipPlane instance exists _dict_properties_factory = { 'slab_separation': lambda self: round(self.slab_separation, 6) + 0, 'repeat_distance': lambda self: round(self.repeat_distance, 6), 'offset': lambda self: round(self.offset, 6) + 0, 'orientation': lambda self: self.orientation.hkl, 'perpendicular_slip_planes': lambda self: list(perp.orientation.hkl for perp in self.perpendicular_slip_planes), 'h_bonded': lambda self: self.h_bonded }
[docs] def as_dict(self): '''Return properties of this slip plane as a dict.''' property_factory = SlipPlanes.SlipPlane._dict_properties_factory return {key: property_factory[key](self) for key in property_factory}
def __init__(self, crystal): '''Initialise a SlipPlanes with a crystal. ''' self._crystal = crystal calculator = ParticleLib.SlipPlaneCalculator(crystal._crystal) self._planes_results = calculator.slip_planes() self.set_separation_threshold(0.0) if no_coordinates(crystal): raise RuntimeError(ParticleLib.basic_critical_text(ParticleLib.no_coordinates_detail_text())) if unknown_spacegroup(crystal): raise RuntimeError(ParticleLib.basic_critical_text(ParticleLib.unknown_spacegroup_detail_text())) check_for_siteless_atoms_or_unknown_bond_types(crystal, warnings, "slip planes") def __getitem__(self, item): return self._planes_list[item] def __len__(self): return len(self._planes_list)
[docs] def separation_threshold(self): '''Return the current separation threshold.''' return self._planes_results.slab_separation_limit()
[docs] def set_separation_threshold(self, minimum_separation): '''Set the separation threshold (in Angstroms), thereby affecting the set of slip planes visible.''' minval = -4.0 if minimum_separation < minval: msg = ( f"Minimum separation threshold {minimum_separation}A is less than minimum limit {minval}A. " f"Actual value used in calculation will be {minval}A." ) warnings.warn(msg) minimum_separation = minval self._planes_results.set_slab_separation_limit(minimum_separation) self._planes_list = [self.SlipPlane(self, plane) for plane in self._planes_results.slip_plane_properties()]
@property def hbond_dimensionality(self): '''Return the hydrogen bond dimensionality as an integer. Returns: -1 no hydrogen bonds, 0 0D discrete, 1 1D chain, 2 2D sheet, 3 3D network ''' return self._planes_results.hbond_dimensionality_int()
[docs] def as_dict(self): '''Return the properties of the slip planes as a dict of lists of values.''' planes = [plane.as_dict() for plane in self] keys = SlipPlanes.SlipPlane._dict_properties_factory.keys() return {key: [plane[key] for plane in planes] for key in keys}