#
# This code is Copyright (C) 2022 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.morphology` module contains classes for crystal morphology.
The main classes in the :mod:`ccdc.morphology` module are:
- :class:`ccdc.morphology.BFDHMorphology`.
- :class:`ccdc.morphology.VisualHabitMorphology`.
'''
import os
import warnings
from dataclasses import dataclass
from enum import IntEnum
from typing import Tuple, Optional, List
from ccdc.molecule import Coordinates, Molecule
from ccdc.crystal import Crystal
from ccdc.utilities import bidirectional_dict, nested_class
from ccdc.utilities import _import_ccdc_modules
_import_ccdc_modules([
'ChemistryLib',
'MathsLib',
'MorphologyLib',
'UtilitiesLib'
])
# Type aliases for crystallographic vectors
IntegerTuple = Tuple[int, int, int]
FloatTuple = Tuple[float, float, float]
#################################################################################################
# Morphology Base class
#################################################################################################
class MorphologyBase:
'''The morphology of a crystal'''
def __init__(self, crystal=None, morph=None):
if crystal is None:
self.identifier = 'DEFAULT'
else:
self.identifier = crystal.identifier
self._morph = morph
self.crystal = crystal
@staticmethod
def from_file(file_name):
'''Creates a Morphology instance from a CIF file.
The CIF file should be those written by this class or Mercury, which
includes a scaling for each of the perpendicular distances.
'''
morph = MorphologyBase()
stream = MorphologyLib.ifstream(file_name)
cif_morph = MorphologyLib.CifMorphology()
morph._morph = cif_morph.read_cif_file(stream)
morph.identifier = os.path.splitext(os.path.basename(file_name))[0]
stream.close()
return morph
@staticmethod
def from_growth_rates(crystal, growth_rates):
'''Creates a morphology from an iterable of growth rates.
:param crystal: an instance of :class:`ccdc.crystal.Crystal`.
:param growth_rates: an iterable of pairs, :class:`ccdc.crystal.Crystal.MillerIndices` and perpendicular distance, otherwise known as
morphological importance.
'''
morph = MorphologyBase()
largest = max(d for mi, d in growth_rates)
morph._morph = MorphologyLib.Morphology(crystal._crystal.cell(), [
MorphologyLib.MorphologicalImportance(mi._miller_indices, d / largest) for mi, d in
growth_rates])
morph.crystal = crystal
morph.identifier = crystal.identifier
return morph
def check_morphology(self):
'''Helper func to check that the morphology has been calculated.'''
if self._morph is None:
raise RuntimeError('No morphology has been calculated.')
def write(self, file_name, keep_all_indices=False):
'''Write this morphology to CIF file.'''
self.check_morphology()
stream = MorphologyLib.ofstream(file_name)
MorphologyLib.CifMorphology().write_cif_file(
stream, self._morph, self.identifier, keep_all_indices)
def relative_area(self, miller_indices):
'''The relative area of the facet.
This is what is usually called the Morphological Importance of a facet.
'''
self.check_morphology()
return self._morph.relative_area(miller_indices._miller_indices)
@property
def facets(self):
'''The facets making up the morphology.'''
self.check_morphology()
facets = self._morph.calculate()
importances = self._morph.morphological_importance_list()
perps = []
for f in facets:
for x in importances:
if x.miller_indices() == f.miller_indices():
perps.append(x.perpendicular_distance() * self._morph.scale_factor())
break
return tuple(
Facet(
f,
p,
Crystal.MillerIndices(
f.miller_indices().h(), f.miller_indices().k(), f.miller_indices().l(),
self.crystal, _cell=self._morph.unit_cell()
)
)
for f, p in zip(facets, perps)
)
@property
def centre_of_geometry(self):
'''The centroid of the morphology.'''
self.check_morphology()
cs = MorphologyLib.SolidAsConvexPolygonSolid(self._morph.solid())
p = cs.centroid()
return Coordinates(p.x(), p.y(), p.z())
@property
def bounding_box(self):
'''The bounding box of the morphology.
A pair of :class:`ccdc.molecule.Coordinates` representing the minimum and maximum corners of the box.
'''
self.check_morphology()
box = self._morph.solid().bounding_box()
return (
Coordinates(box.xmin(), box.ymin(), box.zmin()),
Coordinates(box.xmax(), box.ymax(), box.zmax())
)
@property
def oriented_bounding_box(self):
'''The minimum volume box of the morphology.
This will not necessarily be aligned to the orthonormal cartesian axes.
'''
self.check_morphology()
return OrientedBoundingBox(self)
@property
def volume(self):
'''The volume of the morphology.
This is calculated stochastically, rather than analytically, so has some error.
'''
self.check_morphology()
if not hasattr(self, '_volume'):
volcalc = MorphologyLib.StochasticVolumeCalculator(self._morph.solid())
self._volume = volcalc.solid_volume()
return self._volume
@property
def scale_factor(self):
'''The factor by which the morphology is scaled.'''
self.check_morphology()
return self._morph.scale_factor()
@scale_factor.setter
def scale_factor(self, factor):
'''Set the scale factor of the morphology.'''
self.check_morphology()
self._morph.set_scale_factor(factor)
if hasattr(self, '_volume'):
delattr(self, '_volume')
class OrientedBoundingBox:
'''The bounding box of the morphology.
This box is not necessarily axis-aligned.
'''
def __init__(self, morphology):
self.morphology = morphology
all_pts = [p for f in self.morphology.facets for p in f.coordinates]
points = [MathsLib.Point(c[0], c[1], c[2]) for c in all_pts]
try:
# Fast but may fail with an exception
self._bbox = MathsLib.MinimumVolumeBox(points, MathsLib.MinimumVolumeBox.CALC_REAL)
except RuntimeError:
# Slow and requires lots of stack space
self._bbox = MathsLib.MinimumVolumeBox(points, MathsLib.MinimumVolumeBox.CALC_RATIONAL)
corners = self.corners
from ccdc import descriptors
self.lengths = [
descriptors.GeometricDescriptors.Vector.from_points(corners[0], corners[1]).length,
descriptors.GeometricDescriptors.Vector.from_points(corners[0], corners[2]).length,
descriptors.GeometricDescriptors.Vector.from_points(corners[0], corners[4]).length
]
self.lengths.sort(reverse=True)
@property
def major_length(self):
'''The length of the major axis of the bounding box.'''
return self.lengths[0]
@property
def median_length(self):
'''The length of the middle axis of the bounding box.'''
return self.lengths[1]
@property
def minor_length(self):
'''The length of the minor axis of the bounding box.'''
return self.lengths[2]
@property
def volume(self):
'''The volume of the bounding box.'''
return self._bbox.volume()
@property
def corners(self):
'''The eight points forming the corners of the bounding box.'''
return tuple(
Coordinates(c.x(), c.y(), c.z()) for c in self._bbox.corners()
)
class Facet:
'''One of the facets of a morphology.'''
def __init__(self, _facet, _perpendicular_distance, _miller_indices):
self._facet = _facet
self._miller_indices = _miller_indices
self._perpendicular_distance = _perpendicular_distance
@property
def centre_of_geometry(self):
'''The centre of geometry of the facet.'''
if self._facet is None:
return None
p = self._facet.convex_polygon().center_of_mass_3D()
return Coordinates(p.x(), p.y(), p.z())
@property
def coordinates(self):
'''The coordinates of the vertices of the facet.'''
if self._facet is None:
return None
coords = self._facet.convex_polygon().vertices3D()
return tuple(
Coordinates(*x) for x in coords
)
@property
def edges(self):
'''The edges making up the facet.'''
if self._facet is None:
return None
edges = self._facet.convex_polygon().edges()
return tuple(
(
Coordinates(e.point1().x(), e.point1().y(), e.point1().z()),
Coordinates(e.point2().x(), e.point2().y(), e.point2().z())
)
for e in edges
)
@property
def plane(self):
'''The plane of the facet.
This is a :class:`ccdc.descriptors.GeometricDescriptors.Plane` instance.
'''
if self._facet is None:
return None
from ccdc import descriptors
return descriptors.GeometricDescriptors.Plane(None, None,
_plane=self._facet.convex_polygon().plane())
@property
def miller_indices(self):
'''The Miller indices of the facet.'''
return self._miller_indices
@property
def perpendicular_distance(self):
'''The perpendicular distance from the origin.'''
return self._perpendicular_distance
@property
def area(self):
'''The area of the polygon.'''
if self._facet is None:
return None
return self._facet.convex_polygon().area()
#################################################################################################
# BFDHMorphology
#################################################################################################
[docs]
class BFDHMorphology(MorphologyBase):
_import_ccdc_modules(['BFDHMorphologyLib'])
_telemetry = 0
def __init__(self, crystal=None):
if type(self)._telemetry == 0:
UtilitiesLib.ccdc_bfdh_morphology_telemetry()
type(self)._telemetry = 1
bfdh = BFDHMorphologyLib.BFDHMorphologyCalculator()
morph = bfdh.calculate_morphology(crystal._crystal.cell())
super().__init__(crystal, morph)
#################################################################################################
# VisualHabitMorphology
#################################################################################################
[docs]
class VisualHabitMorphology(MorphologyBase):
'''The morphological information of the VisualHabit calculation'''
def __init__(self, _results, crystal):
self._results = _results
self._morph = self._results.get_morphology()
self.crystal = crystal
self._face_properties = {}
for fp in self._results.get_face_properties():
miller_indices = Crystal.MillerIndices(
fp.miller_indices().h(),
fp.miller_indices().k(),
fp.miller_indices().l(),
self.crystal
)
self._face_properties[miller_indices.hkl] = fp
[docs]
class VisualHabitFacet(Facet):
'''One of the facets of a VisualHabit morphology.'''
def __init__(self, _facet, _perpendicular_distance, _miller_indices, _face_property):
super().__init__(_facet, _perpendicular_distance, _miller_indices)
self._face_property = _face_property
@property
def attachment_energy(self):
'''Attachment energy of the facet in kJ/mol'''
return self._face_property.attachment_energy().get_total()
@property
def surface_energy(self):
'''Surface energy of the facet in mJ/m^2'''
return self._face_property.surface_energy().get_total()
@property
def d_spacing(self):
'''d-spacing of the facet in Angstroms.'''
return self._face_property.d_spacing()
@property
def percentage_area(self):
'''Percentage of the total surface area accounted for by this facet'''
return self._face_property.get_fractional_surface_area() * 100
@property
def offset(self):
'''The offset at which the energies are measured, in Angstroms.'''
return self._face_property.slice_shift()
@property
def facets(self):
'''The facets making up the morphology
:return: A tuple of :class:`ccdc.morphology.VisualHabitMorphology.VisualHabitFacet` instances
'''
facets = self._morph.calculate()
importances = self._morph.morphological_importance_list()
perps = []
for f in facets:
for x in importances:
if x.miller_indices() == f.miller_indices():
perps.append(x.perpendicular_distance() * self._morph.scale_factor())
break
vh_facets = []
for f, p in zip(facets, perps):
miller_indices = Crystal.MillerIndices(
f.miller_indices().h(), f.miller_indices().k(), f.miller_indices().l(),
self.crystal, _cell=self._morph.unit_cell()
)
vh_facets.append(
VisualHabitMorphology.VisualHabitFacet(
f,
p,
miller_indices,
self._face_properties[miller_indices.hkl]
)
)
return tuple(vh_facets)
#################################################################################################
# VisualHabit
#################################################################################################
[docs]
class VisualHabit:
'''
Descriptors for VisualHabit.
'''
_import_ccdc_modules(['HabitMorphologyLib'])
[docs]
class Settings:
'''
Settings for the VisualHabit runner.
'''
_import_ccdc_modules(['HabitMorphologyLib'])
_potentials = bidirectional_dict(
dreidingII=HabitMorphologyLib.VisualHabitPotentialDreiding,
gavezzotti=HabitMorphologyLib.VisualHabitPotentialGavezzotti,
momany=HabitMorphologyLib.VisualHabitPotentialMomany,
CLP=HabitMorphologyLib.VisualHabitPotentialCLP,
CSDOPCS16=HabitMorphologyLib.VisualHabitPotentialCSDOPCS16,
)
def __init__(self):
self._settings = HabitMorphologyLib.VisualHabitRunParameters()
# Default values
_pot = self._potentials.prefix_lookup('dreidingII')
self._settings.set_potential(_pot())
self._settings.set_limiting_radius(30.0)
@property
def potential(self):
'''The potential forcefield used for the calculation.
Set this to one of ``dreidingII`` (default), ``gavezzotti``,
``momany``, ``CLP`` or ``CSD-OPCS16``.
'''
return self._settings.get_potential().information().name()
@potential.setter
def potential(self, value):
pot = self._potentials.prefix_lookup(value.lower().replace('-', ''))
if pot is None:
raise TypeError('Potential was not recognised.')
self._settings.set_potential(pot())
@property
def electrostatic_correction(self):
'''The electrostatic correction mode.
Set to ``None`` to turn this off. Default to ``Evjen``.
'''
return self._settings.get_electrostatic_correction_mode()
@electrostatic_correction.setter
def electrostatic_correction(self, value: float):
self._settings.set_electrostatic_correction_mode(value)
@property
def convergence_limiting_radius(self):
'''The convergence limiting radius for the calculation.
The default is 30.0 Angstroms.
Setting it higher will significantly increase the time of calculation.
'''
return self._settings.get_limiting_radius()
@convergence_limiting_radius.setter
def convergence_limiting_radius(self, value: float):
if value % HabitMorphologyLib.VisualHabitRunParameters.RADIUS_INCREMENT > 0.0:
warnings.warn(
f'The supplied limiting radius of {value} is not a multiple of 2 Å, and will be rounded up.',
UserWarning)
self._settings.set_limiting_radius(float(value))
[docs]
class Results:
'''Holds the results of a VisualHabit calculation.
All energy terms are given in kJ/mol.
'''
_import_ccdc_modules(['HabitMorphologyLib'])
def __init__(self, crystal, _results):
self.crystal = crystal
self._results = _results
self._breakdown = _results.lattice_energies_by_value()
@property
def morphology(self):
'''The calculated morphology
:return: A :class:`ccdc.morphology.VisualHabitMorphology` instance
'''
return VisualHabitMorphology(self._results, self.crystal)
@property
def lattice_energy(self) -> 'VisualHabit.Results.LatticeEnergy':
'''The energies associated with the lattice.
:return: A :class:`ccdc.descriptors.CrystalDescriptors.VisualHabit.Results.LatticeEnergy` instance
'''
return self.LatticeEnergy(
HabitMorphologyLib.VisualHabitEnergyBreakdown(
self._results.total_vdw_attractive_lattice_energy(),
self._results.total_vdw_repulsive_lattice_energy(),
self._results.total_electrostatic_lattice_energy(),
self._results.total_h_bond_attractive_lattice_energy(),
self._results.total_h_bond_repulsive_lattice_energy()
)
)
[docs]
def facet_energy(self, miller_indices: 'Crystal.MillerIndices | tuple'):
'''Return the facet for the miller indices with energy information, but no morphology geometry information.
:param miller_indices: a :class:`ccdc.crystal.Crystal.MillerIndices` instance, or a tuple of (h,k,l).
:return: A :class:`ccdc.morphology.VisualHabitMorphology.VisualHabitFacet` instance.
'''
if isinstance(miller_indices, Crystal.MillerIndices):
mi = miller_indices._miller_indices
else:
mi = ChemistryLib.MillerIndices(*miller_indices)
miller_indices = Crystal.MillerIndices(mi.h(), mi.k(), mi.l(), self.crystal)
props = self._results.get_face_properties(mi)
return VisualHabitMorphology.VisualHabitFacet(None, None, miller_indices, props)
@property
def synthons(self):
'''The synthons, or molecular interactions within the crystal.
:return: A tuple of :class:`ccdc.descriptors.CrystalDescriptors.VisualHabit.Results.Synthon` instances.
'''
synthon_data_store = self._results.synthon_data_store()
synthons = synthon_data_store.get_synthons()
return tuple(
self.Synthon(s, self._results) for s in synthons
)
[docs]
class LatticeEnergy:
'''The lattice energy associated with a facet of the lattice.'''
def __init__(self, _energies):
self._energies = _energies
@property
def total(self) -> float:
'''The total energy.'''
return self._energies.get_total()
@property
def electrostatic(self) -> float:
'''The electrostatic energy.'''
return self._energies.electrostatic()
@property
def h_bond_attraction(self) -> float:
'''The attractive hydrogen bond energy.'''
return self._energies.h_bond_attraction()
@property
def h_bond_repulsion(self) -> float:
'''The repulsive hydrogen bond energy.'''
return self._energies.h_bond_repulsion()
@property
def h_bond(self) -> float:
'''The total hydrogen bond energy.'''
return self.h_bond_attraction + self.h_bond_repulsion
@property
def vdw_attraction(self) -> float:
'''The attractive Van der Waals energy.'''
return self._energies.vdw_attraction()
@property
def vdw_repulsion(self) -> float:
'''The repulsive Van der Waals energy.'''
return self._energies.vdw_repulsion()
@property
def vdw(self) -> float:
'''The total Van der Waals energy.'''
return self.vdw_attraction + self.vdw_repulsion
[docs]
class Synthon:
'''A synthon or molecular interaction within the crystal lattice.'''
def __init__(self, _synthon, _results):
self._synthon = _synthon
self._energies = _synthon.energies()
self._results = _results
self._molecule1 = None
self._molecule2 = None
@property
def distance(self) -> float:
'''The distance between the two molecules's centroids in Angstroms.'''
return self._synthon.separation_ctr_geometry()
@property
def electrostatic(self) -> float:
'''The electrostatic energy in kJ/mol.'''
return self._energies.electrostatic()
@property
def vdw(self) -> float:
'''The Van der Waals energy in kJ/mol.'''
return self._energies.get_vdw_total()
@property
def h_bond(self) -> float:
'''The hydrogen bond energy in kJ/mol.'''
return self._energies.get_h_bonding_total()
@property
def total(self) -> float:
'''The total energy in kJ/mol.'''
return self._energies.get_total()
@property
def symmetry_operator(self) -> str:
'''Return a string representing the symmetry operator relating the two molecules.'''
return self._synthon.distant_transformation().to_string()
@property
def translation_vector(self):
'''Return the unit cell translation vector between the two molecules.'''
offset = self._synthon.distant_cell_offset()
return (offset.x(), offset.y(), offset.z())
@property
def molecule1(self) -> Molecule:
'''Return the first molecule in the synthon.'''
if self._molecule1 is None:
mol = self._results.origin_molecule(self._synthon)
mapping = self._synthon.origin_molecule_index_mapping()
self._molecule1 = Molecule(mapping.label, _molecule=mol)
return self._molecule1
@property
def molecule2(self) -> Molecule:
'''Return the second molecule in the synthon.'''
if self._molecule2 is None:
mol = self._results.distant_molecule(self._synthon)
mapping = self._synthon.distant_molecule_index_mapping()
self._molecule2 = Molecule(mapping.label, _molecule=mol)
return self._molecule2
def __init__(self, settings=None):
if settings is None:
settings = VisualHabit.Settings()
self.settings = settings
[docs]
def calculate(self, crystal: 'Crystal') -> 'VisualHabit.Results':
'''Calculate the habit of the crystal.
:param crystal: a :class:`ccdc.crystal.Crystal` instance.
:returns: a :class:`ccdc.descriptors.CrystalDescriptors.VisualHabit.Results` instance.
The calculation will not be possible if the crystal contains
significant disorder or if some heavy atoms of the crystal have no
coordinates.
'''
csv = crystal._csv
self.settings._settings.set_crystal_view(csv)
if HabitMorphologyLib.check_unknown_bonds(csv):
raise RuntimeError('VisualHabit analysis can not be performed on this structure due to unknown bonds.')
if HabitMorphologyLib.check_disorder_assemblies(csv):
raise RuntimeError(
'VisualHabit analysis can not be performed on this structure due to all or suppressed disorder assemblies selected.')
checker = HabitMorphologyLib.VisualHabitDueDiligence(self.settings._settings)
if not checker.check(csv):
raise RuntimeError(checker.errors_description())
checker.type_atoms()
if not checker.add_charges(self.settings._settings.get_potential()):
raise RuntimeError('VisualHabit analysis can not assign charges to the atoms.')
results = HabitMorphologyLib.VisualHabitRunner(self.settings._settings).run()
return self.Results(crystal, results)
#################################################################################################
# Symmetry Allowed morphology
#################################################################################################
class SymmetryAllowedMorphology(MorphologyBase):
_import_ccdc_modules(['SeedMorphologyLib'])
_telemetry = 0
class Settings:
"""
Settings for the Symmetry Allowed Morphology calculation.
"""
def __init__(self, maximum_hkl_integer: int = 3, minimum_d_spacing: float = 0.5):
self._maximum_hkl_integer = maximum_hkl_integer
self._minimum_d_spacing = minimum_d_spacing
@property
def maximum_hkl_integer(self) -> int:
"""The maximum absolute value of the Miller indices to consider"""
return self._maximum_hkl_integer
@maximum_hkl_integer.setter
def maximum_hkl_integer(self, value: int):
if value < 1:
raise ValueError('maximum_hkl_integer must be a positive integer')
self._maximum_hkl_integer = int(value)
@property
def minimum_d_spacing(self) -> float:
"""The minimum d-spacing to consider"""
return self._minimum_d_spacing
@minimum_d_spacing.setter
def minimum_d_spacing(self, value: float):
if value <= 0.0:
raise ValueError('minimum_d_spacing must be a positive number')
self._minimum_d_spacing = float(value)
def __init__(self, crystal,
settings: 'SymmetryAllowedMorphology.Settings' = None):
"""
Initialize the Symmetry Allowed Morphology object with a crystal structure and settings outlined.
"""
if type(self)._telemetry == 0:
UtilitiesLib.ccdc_seed_morphology_telemetry()
type(self)._telemetry = 1
# TODO: These will be used later when I've implemented the hkl limits
self.settings = settings
if settings is None:
self.settings = self.Settings()
self._crystal = crystal
self._default_growth_rate = 0.5
morph = MorphologyLib.Morphology(crystal._crystal.cell(),
self._default_growth_rate,
self.settings.maximum_hkl_integer)
super().__init__(crystal, morph)
#################################################################################################
# Seed morphology
#################################################################################################
def _vector_to_tuple(vector) -> FloatTuple:
"""Convert a C++ vector to a Python tuple of floats."""
return (vector.x(), vector.y(), vector.z())
def _int_vector_to_tuple(vector) -> IntegerTuple:
"""Convert a C++ integer vector to a Python tuple of integers."""
return (int(vector.x()), int(vector.y()), int(vector.z()))
@dataclass(frozen=True)
class GrowthEnergies:
"""Energy breakdown for a growth element.
All energy values are in kJ/mol.
"""
terrace: float = 0.000
reverse_terrace: float = 0.000
step_edge: float = 0.000
reverse_step_edge: float = 0.000
kink: float = 0.000
reverse_kink: float = 0.000
internal: float = 0.000
bulk: float = 0.000
interfacial: float = 0.000
class SeedMorphology(MorphologyBase):
"""
Morphology based on Hartman-Perdok theory using Periodic Bond Chains (PBCs).
This class models crystal growth using a hierarchical structure:
- **GrowthFace**: A crystal face defined by Miller indices (hkl)
- **Layer**: A repeating growth unit perpendicular to the face
- **StepEdge**: A growth step along a PBC direction
- **EdgeRow**: Individual molecular rows within the step
- **KinkSite**: Attachment sites at row terminations
"""
_import_ccdc_modules(['SeedMorphologyLib'])
_telemetry = 0
@nested_class('SeedMorphology')
class Settings(Crystal.PeriodicBondChain.Settings):
"""Settings for controlling the SeedMorphology calculation."""
def __init__(self):
super().__init__()
@nested_class('SeedMorphology')
class GrowthFace:
"""
A crystal face defined by Miller indices.
Represents a growth face containing layers, step edges, edge rows and kink sites.
"""
@nested_class('GrowthFace')
class GrowthStatusState(IntEnum):
"""Status of a growth calculation element."""
GOOD = 0
BAD = 1
UNFINISHED = 2
def __str__(self) -> str:
return {self.GOOD: "Good", self.BAD: "Bad", self.UNFINISHED: "Unfinished"}[self]
@nested_class('GrowthFace')
class GrowthStatus:
"""Status information for a growth calculation element."""
def __init__(self, status):
"""Initialize the GrowthStatus object."""
self._status = status
def __str__(self) -> str:
message = f"GrowthStatus: Status:{self.status}"
if self.is_good:
return message
return f"{message} Reason:{self.reason})"
@property
def status(self) -> str:
"""Return the status as a human-readable string."""
state = SeedMorphology.GrowthFace.GrowthStatusState(self._status.status())
return str(state)
@property
def reason(self) -> str:
"""Return the reason for a bad status."""
return self._status.reason()
@property
def is_good(self) -> bool:
"""Return True if the status is good."""
return self._status.is_good()
@property
def is_bad(self) -> bool:
"""Return True if the status is bad."""
return self._status.is_bad()
@nested_class('GrowthFace')
class KinkSite:
"""
A kink site on an edge row.
Kink sites are the primary attachment points for molecules during
crystal growth. They represent positions where a molecule can attach
with maximum coordination to the existing crystal surface.
"""
@nested_class("KinkSite")
class _KinkSitePosition(IntEnum):
"""Relative position of a molecule to the kink site."""
AHEAD_OF_KINK = 0
AT_KINK = 1
BEHIND_KINK = 2
NOT_IN_KINK_PBC = 3
def __str__(self) -> str:
return {
self.AHEAD_OF_KINK: "Ahead of kink",
self.AT_KINK: "At kink",
self.BEHIND_KINK: "Behind kink",
self.NOT_IN_KINK_PBC: "Not in kink PBC",
}[self]
class _KinkDirection(IntEnum):
"""Direction of kink propagation."""
EAST = 0
WEST = 1
def __str__(self) -> str:
return {self.EAST: "East", self.WEST: "West"}[self]
def __init__(self, _kink_site):
"""Initialize the KinkSite object."""
self._kink_site = _kink_site
def __str__(self) -> str:
return f"KinkSite: Direction:{self.kink_direction}, Status:{self.growth_status.status}"
@property
def growth_status(self) -> 'GrowthStatus':
"""Return the calculation status for this kink site."""
return SeedMorphology.GrowthFace.GrowthStatus(self._kink_site.status())
def _to_molecule(self, molecule_address) -> Optional[Molecule]:
"""Convert a molecule address to a Molecule object"""
if molecule_address is None:
return None
mol = self._kink_site.interactions_calculation().addressing().obtain_molecule(molecule_address)
return Molecule(_molecule=mol)
@property
def kink_vacancy_position(self) -> FloatTuple:
"""Return the kink vacancy centroid in orthogonal coordinates."""
vacancy = self._kink_site.kink_vacancy()
centroid_point = self._kink_site.interactions_calculation().addressing().molecule_centroid_orth(vacancy)
return _vector_to_tuple(centroid_point)
@property
def preceding_molecule(self) -> Optional[Molecule]:
"""Return the molecule preceding the kink vacancy in the growth direction."""
return self._to_molecule(self._kink_site.preceding_molecule())
@property
def kink_direction_vector(self) -> IntegerTuple:
"""Return the kink growth direction relative to the unit cell. (UVW)"""
return _int_vector_to_tuple(self._kink_site.kink_growth_direction())
@property
def growing_row(self) -> 'Crystal.PeriodicBondChain':
"""Return the PBC of the growing edge row."""
return Crystal.PeriodicBondChain(self._kink_site.growing_row())
@property
def parent_row(self) -> 'EdgeRow':
"""Return the parent edge row of the kink site"""
return SeedMorphology.GrowthFace.EdgeRow(self._kink_site.parent_row())
@property
def offset_kink_x(self) -> float:
"""Return the offset of the kink site in x direction"""
return self._kink_site.offset_kink_x()
def _relative_molecule_position(self, molecule: Molecule) -> str:
"""Return the KinkSitePosition enum as a string for the given molecule address."""
molecule_address = self._kink_site.interactions_calculation().addressing().molecule_address(
molecule._molecule)
pos = self._kink_site.relative_molecule_position(molecule_address)
# SWIG may return an enum proxy or an int; normalise to int then wrap
try:
value = int(pos)
except (TypeError, ValueError):
value = pos
position = SeedMorphology.GrowthFace.KinkSite._KinkSitePosition(value)
return str(position)
@property
def kink_direction(self) -> str:
"""The direction of the kink site with respect ot the parent edge row"""
return str(SeedMorphology.GrowthFace.KinkSite._KinkDirection(self._kink_site.spiral_direction()))
@property
def energies(self) -> GrowthEnergies:
"""
Return all energies for this kink site as a GrowthEnergies object.
> Note - Internal energies are zero for Kink sites
"""
_ks = self._kink_site
return GrowthEnergies(
terrace=_ks.terrace_energy().value(),
reverse_terrace=_ks.reverse_terrace_energy().value(),
step_edge=_ks.step_edge_energy().value(),
reverse_step_edge=_ks.reverse_step_edge_energy().value(),
kink=_ks.kink_growth_energy().value(),
reverse_kink=_ks.reverse_kink_growth_energy().value(),
internal=0.0, # Not available at kink level
bulk=_ks.bulk_energy().value(),
interfacial=_ks.interfacial_energy().value()
)
@nested_class('GrowthFace')
class EdgeRow:
"""
An individual molecular row within a step edge.
Edge rows represent the fundamental growth units that are completed
by kink site attachment during spiral growth.
"""
@nested_class('EdgeRowPosition')
class _EdgeRowPosition(IntEnum):
"""Relative position of a molecule to the edge row."""
BEHIND_CURRENT_ROW = 0
IN_CURRENT_ROW = 1
AHEAD_OF_CURRENT_ROW = 2
def __str__(self) -> str:
return {
self.BEHIND_CURRENT_ROW: "Behind current row",
self.IN_CURRENT_ROW: "In current row",
self.AHEAD_OF_CURRENT_ROW: "Ahead of current row",
}[self]
def __init__(self, _edge_row):
"""Initialize the EdgeRow object."""
self._edge_row = _edge_row
def __str__(self) -> str:
return f"EdgeRow: Edge direction ({self.step_direction}), Status: {self.growth_status.status})"
@property
def growth_status(self) -> 'GrowthStatus':
"""Return the calculation status for this edge row."""
return SeedMorphology.GrowthFace.GrowthStatus(self._edge_row.status())
@property
def parent_step_edge(self) -> 'StepEdge':
"""Return the parent step edge of the edge row - Nice to have"""
return SeedMorphology.GrowthFace.StepEdge(self._edge_row.parent_step_edge())
@property
def step_direction(self) -> IntegerTuple:
"""Return the step direction of the edge row."""
return _int_vector_to_tuple(self._edge_row.step_direction())
@property
def _current_row(self) -> 'Crystal.PeriodicBondChain':
"""Return the current row of the edge row"""
return Crystal.PeriodicBondChain(self._edge_row.current_row(),
self._edge_row.interactions_calculation())
@property
def _next_row(self) -> 'Crystal.PeriodicBondChain':
"""Return the next row of the edge row"""
return Crystal.PeriodicBondChain(self._edge_row.next_row(), self._edge_row.interactions_calculation())
@property
def parent_face(self) -> 'GrowthFace':
"""Return the parent growth face of the edge row - Nice to have"""
return SeedMorphology.GrowthFace(self._edge_row.parent_face())
@property
def row_orth_normal(self) -> FloatTuple:
"""Return the orthogonal row normal of the edge row."""
return _vector_to_tuple(self._edge_row.row_orth_normal())
@property
def row_frac_normal(self) -> FloatTuple:
"""Return the fractional row normal of the edge row."""
return _vector_to_tuple(self._edge_row.row_frac_normal())
@property
def row_orth_vector(self) -> FloatTuple:
"""Return the vector representing one full growth cycle along the edge row in orthogonal coordinates."""
return _vector_to_tuple(self._edge_row.row_orth_vector())
@property
def average_width_per_growth_unit(self) -> float:
"""
Return the average width per growth unit along the edge row.
"""
return self._edge_row.average_width_per_growth_unit()
@property
def _current_min_step_y(self) -> float:
"""Return the minimum y value of the current edge row"""
return self._edge_row.current_min_step_y()
@property
def _current_max_step_y(self) -> float:
"""Return the maximum y value of the current edge row"""
return self._edge_row.current_max_step_y()
@property
def _centroid_step_y(self) -> float:
"""Return the centroid y value of the edge row"""
return self._edge_row.centroid_step_y()
@property
def kink_sites(self) -> Tuple['KinkSite', ...]:
"""Return the kink sites of the edge row"""
return tuple(
SeedMorphology.GrowthFace.KinkSite(ks)
for ks in self._edge_row.kink_sites()
)
def _next_kink_site(self, kink_site) -> 'KinkSite':
"""Return the next kink site of the edge row"""
ks = self._edge_row.next_kink_site(kink_site._kink_site)
if ks is None:
return None
return SeedMorphology.GrowthFace.KinkSite(ks)
def _previous_kink_site(self, kink_site) -> 'KinkSite':
"""Return the previous kink site of the edge row"""
ks = self._edge_row.previous_kink_site(kink_site._kink_site)
if ks is None:
return None
return SeedMorphology.GrowthFace.KinkSite(ks)
@property
def kink_sites_by_direction(self) -> List[List['KinkSite']]:
"""Return a dictionary of kink sites by growth direction"""
k_sites_by_direction = []
for kink_sites in self._edge_row.kink_sites_by_direction():
k_sites = []
for ks in kink_sites:
k_sites.append(SeedMorphology.GrowthFace.KinkSite(ks))
k_sites_by_direction.append(k_sites)
return k_sites_by_direction
def _relative_molecule_position(self, molecule: Molecule) -> str:
"""Return the relative position of molecules to the edge row."""
molecule_address = self._edge_row.interactions_calculation().addressing().molecule_address(
molecule._molecule)
position = SeedMorphology.GrowthFace.EdgeRow._EdgeRowPosition(
self._edge_row.relative_molecule_position(molecule_address))
return str(position)
@property
def energies(self) -> GrowthEnergies:
"""
Return all energies for this edge row as a GrowthEnergies object.
> Note - Kink and Reverse Kink energies are zero for edge rows
"""
_er = self._edge_row
return GrowthEnergies(
terrace=_er.terrace_energy().value(),
reverse_terrace=_er.reverse_terrace_energy().value(),
step_edge=_er.step_edge_energy().value(),
reverse_step_edge=_er.reverse_step_edge_energy().value(),
internal=_er.internal_energy().value(),
bulk=_er.bulk_energy().value(),
interfacial=_er.interfacial_energy().value(),
)
@nested_class('GrowthFace')
class StepEdge:
"""
A growth step along a PBC direction.
Step edges represent an incomplete layer.
They contain edge rows which in turn contain kink sites.
"""
@nested_class('StepEdgePosition')
class _StepEdgePosition(IntEnum):
"""The relative position of a molecule with respect to the step edge"""
BEHIND_STEP = 0
IN_STEP = 1
AHEAD_OF_STEP = 2
NOT_IN_STEP = 3
def __str__(self) -> str:
return {
self.BEHIND_STEP: "Behind step",
self.IN_STEP: "In step",
self.AHEAD_OF_STEP: "Ahead of step",
self.NOT_IN_STEP: "Not in step",
}[self]
def __init__(self, _step_edge):
"""Initialize the StepEdge object."""
self._step_edge = _step_edge
def __str__(self) -> str:
return f"StepEdge: Step Edge Direction {self.step_direction}, Status: {self.growth_status.status})"
@property
def growth_status(self) -> 'GrowthStatus':
"""Return the calculation status for this step edge."""
return SeedMorphology.GrowthFace.GrowthStatus(self._step_edge.status())
@property
def parent_layer(self) -> 'Layer':
"""Return the parent edge row of the edge row - Nice to have"""
return SeedMorphology.GrowthFace.Layer(self._step_edge.parent_layer())
@property
def step_direction(self) -> IntegerTuple:
"""Return the step direction of the step edge."""
return _int_vector_to_tuple(self._step_edge.step_direction())
@property
def edge_rows(self) -> Tuple['EdgeRow', ...]:
"""Return the edge rows of the step edge."""
return tuple(SeedMorphology.GrowthFace.EdgeRow(er) for er in self._step_edge.edge_rows())
@property
def step_orth_normal(self) -> FloatTuple:
"""Return the orthogonal step normal of the step edge."""
return _vector_to_tuple(self._step_edge.step_orth_normal())
@property
def step_frac_normal(self) -> FloatTuple:
"""Return the fractional step normal of the step edge."""
return _vector_to_tuple(self._step_edge.step_frac_normal())
@property
def average_width_per_growth_unit(self) -> float:
"""
Return the average width per growth unit along the rows.
"""
return self._step_edge.average_width_per_growth_unit()
@property
def average_propagation_length(self) -> float:
"""
Return the average propagation length per row in the step edge direction.
"""
return self._step_edge.average_propagation_length()
@property
def energies(self) -> GrowthEnergies:
"""
Return all energies for this step edge as a GrowthEnergies object.
> Note - Kink and Reverse Kink energies are zero for Step edges"""
_se = self._step_edge
return GrowthEnergies(
terrace=_se.terrace_energy().value(),
reverse_terrace=_se.reverse_terrace_energy().value(),
step_edge=_se.step_edge_energy().value(),
reverse_step_edge=_se.reverse_step_edge_energy().value(),
internal=_se.internal_energy().value(),
bulk=_se.bulk_energy().value(),
interfacial=_se.interfacial_energy().value(),
)
@nested_class('GrowthFace')
class Layer:
"""
A single growth layer within a crystal face.
Represents the fundamental repeating unit perpendicular to the face normal.
Contains step edges that define the layer's growth directions.
In the Hartman-Perdok theory, layers correspond to F-faces (flat faces)
when they contain two or more non-parallel PBCs.
"""
def __init__(self, _layer):
"""Initialize the Layer object."""
self._layer = _layer
if self._layer is None:
raise ValueError('layer cannot be None')
def __str__(self) -> str:
return f"Layer: Status: {self.growth_status.status}"
@property
def growth_status(self) -> 'GrowthStatus':
"""Return the calculation status for this layer."""
return SeedMorphology.GrowthFace.GrowthStatus(self._layer.status())
@property
def parent_face(self) -> 'GrowthFace':
"""Return the parent growth face of the layer"""
return SeedMorphology.GrowthFace(self._layer.parent_face())
@property
def periodic_bond_chains(self) -> List['Crystal.PeriodicBondChain']:
"""Return the periodic_bond_chains of the growth face"""
return [Crystal.PeriodicBondChain(pbc, pbc.interactions_calculation())
for pbc in self._layer.pbcs()]
@property
def molecules(self) -> List[Molecule]:
"""Return the molecule addresses in the layer"""
return [Molecule(_molecule=self._layer.interactions_calculation().addressing().obtain_molecule(addr))
for addr in self._layer.molecule_addresses()]
@property
def min_z_position(self) -> float:
"""Return the minimum z value of the layer - where z is perpendicular to the layer"""
return self._layer.min_face_z()
@property
def max_z_position(self) -> float:
"""Return the maximum z value of the layer - where z is perpendicular to the layer"""
return self._layer.max_face_z()
@property
def centroid_z_position(self) -> float:
"""Return the centroid of the layer - relative to the layer perpendicular"""
return self._layer.centroid_face_z()
@property
def step_edges(self) -> Tuple['StepEdge', ...]:
"""Return the step edges of the layer"""
return tuple(SeedMorphology.GrowthFace.StepEdge(er) for er in self._layer.step_edges())
@property
def step_edge_directions(self) -> List[IntegerTuple]:
"""Return the step edge directions of the layer"""
return [_int_vector_to_tuple(sed)
for sed in self._layer.step_edge_directions()]
def step_edge_by_direction(self, direction: Tuple[int, int, int]) -> 'StepEdge':
"""Return the step edge for the given direction"""
try:
se = self._layer.step_edge_by_direction(MathsLib.int_vector_3d(*direction))
return SeedMorphology.GrowthFace.StepEdge(se)
except RuntimeError as e:
if 'map' in str(e):
raise KeyError(f'Step edge with direction {direction} not found.')
raise
@property
def energies(self) -> GrowthEnergies:
"""
Return all energies for this layer as a GrowthEnergies object.
Note - Step Edge, Reverse Step edge, Kink, Reverse Kink energies are zero for edge rows.
"""
return GrowthEnergies(
terrace=self._layer.terrace_energy().value(),
reverse_terrace=self._layer.reverse_terrace_energy().value(),
internal=self._layer.internal_energy().value(),
bulk=self._layer.reverse_terrace_energy().value(),
interfacial=self._layer.terrace_energy().value(),
)
def __init__(self, growth_face):
"""Initialize the GrowthFace object."""
self._growth_face = growth_face
def __str__(self) -> str:
return f"GrowthFace: MillerIndices{self.miller_indices}, Status={self.status.status}"
@property
def status(self) -> GrowthStatus:
"""Return the calculation status for this growth face."""
return SeedMorphology.GrowthFace.GrowthStatus(self._growth_face.status())
@property
def layers(self) -> Tuple['Layer', ...]:
"""Return the layers of the growth face."""
return tuple(
SeedMorphology.GrowthFace.Layer(layer)
for layer in self._growth_face.layers()
)
@property
def periodic_bond_chains(self) -> Tuple[Crystal.PeriodicBondChain, ...]:
"""Return the periodic bond chains of the growth face."""
return tuple(
Crystal.PeriodicBondChain(pbc, pbc.interactions_calculation())
for pbc in self._growth_face.pbcs()
)
@property
def face_orth_normal(self) -> FloatTuple:
"""Return the orthogonal face normal of the growth face."""
return _vector_to_tuple(self._growth_face.face_orth_normal())
@property
def face_frac_normal(self) -> FloatTuple:
"""Return the fractional face normal of the growth face."""
return _vector_to_tuple(self._growth_face.face_frac_normal())
@property
def num_symmetry_layers(self) -> int:
"""Return the number of layers implied by space group symmetry"""
return self._growth_face.n_symmetry_layers()
@property
def interplanar_spacing(self) -> float:
"""Return the interplanar spacing of the growth face in Angstroms."""
return self._growth_face.interplanar_spacing()
@property
def average_height_per_layer(self) -> float:
"""
Return the average height per layer in Angstroms.
In practice this may differ from the interplanar spacing.
"""
return self._growth_face.average_height_per_layer()
def get_layer_containing_molecule(self, mol: Molecule) -> Optional['Layer']:
"""Return the layer containing the given molecule address"""
mol_address = self._growth_face.interactions().addressing().molecule_address(mol._molecule)
layer = self._growth_face.get_layer_containing_address(mol_address)
if layer is None:
return None
return SeedMorphology.GrowthFace.Layer(layer)
@property
def miller_indices(self) -> IntegerTuple:
"""Return the Miller indices of the growth face."""
return _int_vector_to_tuple(self._growth_face.miller_indices().to_vector())
@nested_class('SeedMorphology')
class GrowthFacet(Facet):
"""
A facet on a SeedMorphology with associated growth face information.
Extends the base Facet class to include a reference to the GrowthFace
which contains the detailed growth calculation results.
"""
def __init__(self, facet, perpendicular_distance, miller_indices, growth_face=None):
"""Initialize the GrowthFacet object with a morphology.Facet object."""
super().__init__(facet, perpendicular_distance, miller_indices)
self._growth_face = growth_face
def __str__(self) -> str:
return f"GrowthFacet: MillerIndices:{self.miller_indices.hkl}"
@property
def growth_face(self) -> Optional['SeedMorphology.GrowthFace']:
"""Return the growth face associated with this facet."""
if self._growth_face is None:
return None
return SeedMorphology.GrowthFace(self._growth_face)
def __init__(self, crystal, settings=None):
"""
Initialize the SeedMorphology object with a crystal structure and settings.
"""
if type(self)._telemetry == 0:
UtilitiesLib.ccdc_seed_morphology_telemetry()
type(self)._telemetry = 1
self._crystal = crystal
self._settings = settings
self._relative_growth_rate = 0.5
if self._settings is None:
self._settings = SeedMorphology.Settings()
morph = None
pbcs = self._crystal.periodic_bond_chains(self._settings)
if len(pbcs) == 0:
warnings.warn(
f"No periodic bond chains found for crystal '{crystal.identifier}'. "
"Morphology properties will be unavailable.",
UserWarning
)
self._growth_faces = []
self._interaction_calculation = None
else:
self._interaction_calculation = pbcs[0]._interaction_calculation
chains = [pbc._pbc for pbc in pbcs]
calculator = SeedMorphologyLib.GrowthFaceCalculator(self._interaction_calculation, chains)
self._growth_faces = calculator.calculate_growth_faces(None)
morph = calculator.calculate_morphology(self._crystal._crystal.cell(), self._relative_growth_rate)
super().__init__(crystal, morph)
@property
def facets(self):
"""
Overloaded function from the MorphologyBase - Return the facets of the seed -
"""
facets = self._morph.calculate()
growth_rate_scaled = self._relative_growth_rate * self._morph.scale_factor()
growth_faces_map = {
(gf.miller_indices().h(), gf.miller_indices().k(), gf.miller_indices().l()): gf
for gf in self._growth_faces
}
growth_facets = []
for facet in facets:
hkl = (facet.miller_indices().h(), facet.miller_indices().k(), facet.miller_indices().l())
growth_facet = SeedMorphology.GrowthFacet(
facet=facet,
perpendicular_distance=growth_rate_scaled,
miller_indices=Crystal.MillerIndices(*hkl, self.crystal, _cell=self._morph.unit_cell()),
growth_face=growth_faces_map.get(hkl, None))
growth_facets.append(growth_facet)
return tuple(growth_facets)
def get_facet(self, hkl: IntegerTuple) -> Optional['GrowthFacet']:
"""Return the facet for the given Miller indices.
:param hkl: Miller indices as (h, k, l) tuple.
:returns: The matching GrowthFacet, or None if not found.
"""
for facet in self.facets:
if facet.miller_indices.hkl == hkl:
return facet
return None
@property
def growth_faces(self) -> Tuple['GrowthFace', ...]:
"""Return all calculated growth faces."""
return tuple(SeedMorphology.GrowthFace(gf) for gf in self._growth_faces)
def get_growth_face(self, hkl: IntegerTuple) -> Optional['GrowthFace']:
"""Return the growth face for the given Miller indices.
:param hkl: Miller indices as (h, k, l) tuple.
:returns: The matching GrowthFace, or None if not found.
"""
for gf in self._growth_faces:
indices = gf.miller_indices()
if (indices.h(), indices.k(), indices.l()) == hkl:
return SeedMorphology.GrowthFace(gf)
return None