#
# 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])
(23.0, -0.0, -0.0)
>>> round(len(list(particle_surface.topology.triangles)), -3)
8000
>>> tuple(round(x,0) for x in list(particle_surface.topology.triangles)[0][0])
(23.0, 0.0, -0.0)
>>> particle_surface.descriptors.projected_area
165.167
>>> particle_surface.descriptors.rugosity
1.375
>>> particle_surface.descriptors.rmsd
1.642
>>> particle_surface.descriptors.skewness
0.777
>>> particle_surface.descriptors.kurtosis
2.793
>>> 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.crystal import Crystal
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`
:param hbond_criterion: A hydrogen bond criterion to use in the surface analysis, defaults to the Mercury default criterion
:type hbond_criterion: :class:`ccdc.molecule.Molecule.HBondCriterion`, optional
'''
[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, hbc):
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, hbc)
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, hbc):
self._a, self._d, self._ud = ParticleLib.get_atom_property_list(cssv, contacts, hbc)
@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,
hbond_criterion=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
if hbond_criterion:
hbc = hbond_criterion._hbc.hbond_criterion()
else:
hbc = molecule.Molecule.HBondCriterion()._hbc.hbond_criterion()
hbc.set_allow_non_hydrogen_contacts(False)
self._descriptors = Surface.Descriptors(_topology,
csv,
self._cssv,
self._contacts,
hbc)
self.atom_properties = Surface.AtomProperties(self._cssv, self._contacts, hbc)
@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`
:param hbond_criterion: A hydrogen bond criterion to use in the surface analysis, defaults to the Mercury default criterion
:type hbond_criterion: :class:`ccdc.molecule.Molecule.HBondCriterion`, optional
'''
[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, hbond_criterion=None):
'''Initialise a SlipPlanes with a crystal and optional hydrogen bond criterion.
'''
if hbond_criterion:
hbc = hbond_criterion._hbc.hbond_criterion()
else:
hbc = molecule.Molecule.HBondCriterion()._hbc.hbond_criterion()
hbc.set_allow_non_hydrogen_contacts(False)
self._crystal = crystal
self._calculator = ParticleLib.SlipPlaneCalculator(crystal._csv, hbc)
self._planes_results = self._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}
[docs] def slip_planes_with_orientation(self, miller_indices):
'''Return all slip planes aligned with the given Miller plane, ordered by greatest separation first.'''
if isinstance(miller_indices, Crystal.MillerIndices):
mi = miller_indices._miller_indices
else:
mi = ChemistryLib.MillerIndices(*miller_indices)
planes = self._calculator.all_planes_on_axis(mi)
return [self.SlipPlane(self, plane) for plane in planes]