#
# 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 ccdc.molecule import Coordinates
from ccdc.crystal import Crystal
from ccdc.utilities import bidirectional_dict, Logger
from ccdc.utilities import _private_importer
with _private_importer() as pi:
pi.import_ccdc_module('ChemistryLib')
pi.import_ccdc_module('MathsLib')
pi.import_ccdc_module('BFDHMorphologyLib')
pi.import_ccdc_module('UtilitiesLib')
#################################################################################################
# Morphology Base class
#################################################################################################
class MorphologyBase:
'''The morphology of a crystal'''
def __init__(self, crystal=None):
if crystal is None:
self._morph = BFDHMorphologyLib.Morphology()
self.identifier = 'DEFAULT'
else:
bfdh = BFDHMorphologyLib.BFDHMorphologyCalculator()
self._morph = bfdh.calculate_morphology(crystal._crystal.cell())
self.identifier = crystal.identifier
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 = BFDHMorphologyLib.ifstream(file_name)
cif_morph = BFDHMorphologyLib.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 = BFDHMorphologyLib.Morphology(crystal._crystal.cell(), [
BFDHMorphologyLib.MorphologicalImportance(mi._miller_indices, d / largest) for mi, d in
growth_rates])
morph.crystal = crystal
morph.identifier = crystal.identifier
return morph
def write(self, file_name, keep_all_indices=False):
'''Write this morphology to CIF file.'''
stream = BFDHMorphologyLib.ofstream(file_name)
BFDHMorphologyLib.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.
'''
return self._morph.relative_area(miller_indices._miller_indices)
@property
def facets(self):
'''The facets making up the 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.'''
cs = BFDHMorphologyLib.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.
'''
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.
'''
return OrientedBoundingBox(self)
@property
def volume(self):
'''The volume of the morphology.
This is calculated stochastically, rather than analytically, so has some error.
'''
if not hasattr(self, '_volume'):
volcalc = BFDHMorphologyLib.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.'''
return self._morph.scale_factor()
@scale_factor.setter
def scale_factor(self, factor):
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.'''
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.'''
coords = self._facet.convex_polygon().vertices3D()
return tuple(
Coordinates(*x) for x in coords
)
@property
def edges(self):
'''The edges making up the facet.'''
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.
'''
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.'''
return self._facet.convex_polygon().area()
#################################################################################################
# BFDHMorphology
#################################################################################################
[docs]class BFDHMorphology(MorphologyBase):
_telemetry = 0
def __init__(self, crystal=None):
super().__init__(crystal)
if type(self)._telemetry == 0:
UtilitiesLib.ccdc_bfdh_morphology_telemetry()
type(self)._telemetry = 1
#################################################################################################
# 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'''
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 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.
'''
with _private_importer() as pi:
pi.import_ccdc_module('HabitMorphologyLib')
[docs] class Settings:
'''
Settings for the VisualHabit runner.
'''
with _private_importer() as pi:
pi.import_ccdc_module('HabitMorphologyLib')
_potentials = bidirectional_dict(
dreidingII = HabitMorphologyLib.VisualHabitPotentialDreiding,
gavezzotti = HabitMorphologyLib.VisualHabitPotentialGavezzotti,
momany = HabitMorphologyLib.VisualHabitPotentialMomany,
)
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``,
or ``momany``.
'''
return self._settings.get_potential().information_.name_
@potential.setter
def potential(self, value):
pot = self._potentials.prefix_lookup(value.lower())
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):
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):
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.
'''
with _private_importer() as pi:
pi.import_ccdc_module('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):
'''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] class LatticeEnergy:
'''The lattice energy associated with a facet of the lattice.'''
def __init__(self, _energies):
self._energies = _energies
@property
def total(self):
'''The total energy.'''
return self._energies.get_total()
@property
def electrostatic(self):
'''The electrostatic energy.'''
return self._energies.electrostatic()
@property
def h_bond_attraction(self):
'''The attractive hydrogen bond energy.'''
return self._energies.h_bond_attraction()
@property
def h_bond_repulsion(self):
'''The repulsive hydrogen bond energy.'''
return self._energies.h_bond_repulsion()
@property
def h_bond(self):
'''The total hydrogen bond energy.'''
return self.h_bond_attraction + self.h_bond_repulsion
@property
def vdw_attraction(self):
'''The attractive Van der Waals energy.'''
return self._energies.vdw_attraction()
@property
def vdw_repulsion(self):
'''The repulsive Van der Waals energy.'''
return self._energies.vdw_repulsion()
@property
def vdw(self):
'''The total Van der Waals energy.'''
return self.vdw_attraction + self.vdw_repulsion
def __init__(self, settings=None):
if settings is None:
settings = VisualHabit.Settings()
self.settings = settings
[docs] def calculate(self, crystal):
'''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 = ChemistryLib.CrystalStructureView_instantiate(crystal._crystal)
self.settings._settings.set_crystal_view(csv)
if HabitMorphologyLib.check_unknown_bonds(crystal._crystal):
raise RuntimeError('VisualHabit analysis can not be performed on this structure due to unknown bonds.')
if HabitMorphologyLib.check_disorder_assemblies(crystal._crystal):
raise RuntimeError('VisualHabit analysis can not be performed on this structure due to disorder assemblies.')
checker = HabitMorphologyLib.VisualHabitDueDiligence(self.settings._settings)
if not checker.check(crystal._crystal):
raise RuntimeError(checker.errors_description())
checker.gasteiger_charges()
checker.type_atoms()
results = HabitMorphologyLib.VisualHabitRunner(self.settings._settings).run()
return self.Results(crystal, results)
