Morphology

There is a class, ccdc.morphology.BFDHMorphology, representing the morphology/habit of a crystal.

To make use of the Morphology class, first we need to import BFDHMorphology, and since we will be visualising the predicted morphology using matplotlib we will import its python interface:

>>> import matplotlib
>>> matplotlib.use('Agg')
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D
>>> from mpl_toolkits.mplot3d.art3d import Poly3DCollection
>>> from ccdc.morphology import BFDHMorphology

In the example script we will look at the morphology of the beta polymorph of succinic-acid: CSD refcode SUCACB03

SUCACB03 2D diagram.

To create a Morphology object of SUCACB03 we read the crystal in the usual way:

>>> from ccdc.io import EntryReader
>>> reader = EntryReader('csd')
>>> crystal = reader.crystal('SUCACB03')

Then we create a morphology object like so:

>>> morphology = BFDHMorphology(crystal)

By default the morphology produced in this way will be a calculated BFDH morphology. BFDH (Bravais, Friedel, Donnay and Harker) morphology is a theoretical morphological prediction based purely on the geometry of the crystal lattice, and is based on the assumption that crystal face growth is inversely proportional to the inter-planar distance. See the article “American Mineralogist (1937), 22, 446-67” for more details on the calculation.

Later we will also look at reading custom, and experimental morphologies into the API

A morphology instance has various attributes that we have access to:

>>> print(morphology.identifier)
SUCACB03
>>> print(morphology.centre_of_geometry) 
Coordinates(x=0.000, y=-0.000, z=0.000)
>>> print(morphology.bounding_box) 
(Coordinates(x=-7.044, y=-8.131, z=-9.451), Coordinates(x=7.044, y=8.131, z=9.451))
>>> print(round(morphology.volume, 4)) 
    2786.5006

The volume of the morphology is calculated by a stochastic method, so contains an amount of error. If this is deemed too great, the sampling frequency can be increased, though there is no API method to do this.

As well as the coordinates of the orthonormal bounding box there is an attribute of the morphology to calculate an oriented bounding box. This will be very close to the minimum enclosing volume of the morphology, and may be used to calculate aspect ratios of the morphology:

>>> bbox = morphology.oriented_bounding_box
>>> print(round(bbox.volume, 3))
4162.907
>>> print(round(bbox.major_length, 3))
19.59
>>> print(round(bbox.median_length, 3))
16.263
>>> print(round(bbox.minor_length, 3))
13.067
>>> print(bbox.corners) 
(Coordinates(x=-5.960, y=8.131, z=-10.154), Coordinates(x=-7.085, y=8.131, z=9.404), Coordinates(x=-5.960, y=-8.131, z=-10.154), Coordinates(x=-7.085, y=-8.131, z=9.404), Coordinates(x=7.085, y=8.131, z=-9.404), Coordinates(x=5.960, y=8.131, z=10.154), Coordinates(x=7.085, y=-8.131, z=-9.404), Coordinates(x=5.960, y=-8.131, z=10.154))

The morphology can also be subdivided into faces, edges and vertices. The faces of the morphology, are instances of the class ccdc.morphology.Facet.

>>> facets = morphology.facets
>>> f = facets[0]
>>> print(f.centre_of_geometry) 
Coordinates(x=6.523, y=-0.000, z=0.375)
>>> print(round(f.perpendicular_distance, 3))
18.332
>>> print(f.coordinates[0]) 
Coordinates(x=6.911, y=4.066, z=-6.384)
>>> print(f.edges[0]) 
(Coordinates(x=6.134, y=4.066, z=7.134), Coordinates(x=6.911, y=4.066, z=-6.384))
>>> print(round(f.area, 3))
128.972
>>> print(f.miller_indices) 
MillerIndices(1, 0, 0)
>>> f.plane 
Plane(normal=Vector(0.998, 0.000, 0.057), distance=6.533)

The complete set of miller indices appropriate to the morphology can be retrieved by:

>>> all_miller_indices = [f.miller_indices for f in facets]
>>> print(all_miller_indices) 
[MillerIndices(1, 0, 0), MillerIndices(-1, 0, 0), MillerIndices(-1, -1, 0), MillerIndices(1, 1, 0), MillerIndices(1, -1, 0), MillerIndices(-1, 1, 0), MillerIndices(0, -2, 0), MillerIndices(0, 2, 0), MillerIndices(0, -1, 1), MillerIndices(0, 1, 1), MillerIndices(0, 1, -1), MillerIndices(0, -1, -1), MillerIndices(-1, 1, 1), MillerIndices(1, 1, -1), MillerIndices(-1, -1, 1), MillerIndices(1, -1, -1)]

We can get the relative area or morphological importance associated with an instance of the ccdc.crystal.Crystal.MillerIndices class:

>>> print(round(morphology.relative_area(all_miller_indices[0]), 3))
0.119

They can be associated with their relative areas:

>>> [round(morphology.relative_area(mi), 3) for mi in all_miller_indices]
[0.119, 0.119, 0.052, 0.052, 0.052, 0.052, 0.071, 0.071, 0.098, 0.098, 0.098, 0.098, 0.005, 0.005, 0.005, 0.005]

The morphology can also be visualised using Matplotlib, and the following visualisation functions can be found in the example script

First, we set up the Matplotlib figure:

>>> fig = plt.figure(figsize=(10., 10.))   # 3D graph instance
>>> ax = fig.add_subplot(111, projection='3d')   # 3D Axes
>>> _ = ax.set_xlim3d(-15, 15)  # Set the axes limits
>>> _ = ax.set_ylim3d(-15, 15)
>>> _ = ax.set_zlim3d(-15, 15)
>>> _ = ax.grid(False)  # To hide the gridlines
>>> _ = plt.axis('off')  # To hide the axes

The following function will plot the edges of each facet given to it into the 3D plot, resulting in a wireframe representation:

>>> def wireframe(facet):
...     for edge in facet.edges:
...         Axes3D.plot(ax,
...                                         [coord[0] for coord in edge],  # The x coordinates
...                     [coord[1] for coord in edge],  # The y coordinates
...                     [coord[2] for coord in edge],  # The z coordinates
...                                 c='black',
...                                         linewidth=1.5)

In order to mesh the wireframe, the following function can be used. Mesh surfaces are made up of triangles, and so the mesh for each crystal face is formed by taking the two coordinates of each edge,along with “centre_of_geometry” of the face to produce the triangles needed:

>>> def skin(facet):
...     for edge in facet.edges:
...         vertices = [(edge[0], edge[1], facet.centre_of_geometry)]
...         Axes3D.add_collection3d(ax, Poly3DCollection(vertices, color='blue', linewidth=0, alpha=0.3))

Finally, labels can be added to the plot, In this case the hkl, and relative area, are plotted at each face centre_of_geometry along with a point (using scatter)

>>> def add_labels(morphology, face):
...     ax.text(face.centre_of_geometry[0], face.centre_of_geometry[1], face.centre_of_geometry[2],
...             '''   {}
...     {}'''.format(face.miller_indices.hkl, round(morphology.relative_area(face.miller_indices), 3)),
...             color='black')
...
...     ax.scatter(face.centre_of_geometry[0],
...                face.centre_of_geometry[1],
...                face.centre_of_geometry[2],
...                s=10,
...                color='black')

To generate the final plot, we give the functions the facets of the morphology (this is pre-written in the function generate_plot in the example):

>>> for facet in morphology.facets:
...     wireframe(facet)
...     skin(facet)
...     add_labels(morphology, facet)

This gives the image:

Morphology of SUCACBO3

Further analysis can be performed, for instance we can calculate the angles between each facet in the morphology using the following function:

>>> def compare_angles(morphology):
...     indices = {}
...     for f in morphology.facets:
...         indices[f.miller_indices.hkl] = GeometricDescriptors.Plane.from_points(*f.coordinates)
...     combinations = list(itertools.combinations([k for k in indices], 2))  # Every pair of facets to compare
...     for combination in combinations:
...         angle = 180 - (indices[combination[0]]).plane_angle(indices[combination[1]])
...         print('planes {} & {}: angle = {}'.format(combination[0], combination[1], angle))

These angles could then be used to interpret optical goniometry data (CrystEngComm, 2015, 17, 1015)

Also included is a function to present a nicely formatted table of relative face areas

>>> def print_areas(morphology):
...     for face in morphology.facets:
...         print('{}: relative area = {}'.format(
...             face.miller_indices.hkl, round(morphology.relative_area(face.miller_indices), 3)))

As a ‘run_all’ function has been produced, the example script will perform all of this automatically.

As mentioned previously BFDH is a simplistic approximation, and so reading morphologies in from more advanced predictions, or experimental observation may be required. To do this, morphologies can be read from morphology cifs (prediction outputs) or CIF files that contain experimentally indexed faces

Supplied with this document is the file ‘experimental.cif’. This is a real experimental structure, where the faces of the crystal have been measured and recorded.

If we change the line:

>>> morphology = BFDHMorphology(crystal)

to use this file:

>>> cif_file = 'experimental.cif'

>>> morphology = BFDHMorphology.from_file(cif_file)

We will instead read in the experimental faces to our morphology object. If the CIF is not located in the same folder as your script, you will need to supply the path as well as the file name. Morphology CIF files such as those save from Mercury’s BFDH calculator can also be read in the same way

Morphology instances can also be written to morphology_cif files with the following lines

>>> import os, tempfile
>>> temp_dir = tempfile.mkdtemp()
>>> file_name = os.path.join(temp_dir, 'AABHTZ_morph.cif')
>>> morphology.write(file_name)

In addition one has the ability to construct a BFDH morphology from an iterable of morphological importances, here provided as pairs of ccdc.crystal.Crystal.MillerIndices and a perpendicular distance. In general the reliable determination of growth rates is stll an unsolved problem due to the complexity of crystallisation mechanisms, but this method allows researchers to experiment with any methods they may have developed, and to visualise the results. Care must be taken to ensure that the generated morphology will subtend a closed convex polyhedron: if not various methods of the morphology will fail.

Here we exemplify its use with some calculated growth rates:

>>> from ccdc.io import EntryReader
>>> csd = EntryReader('csd')
>>> hxacan = csd.crystal('HXACAN')
>>> growth_rates = [
...     (hxacan.MillerIndices(2, 1, 1, hxacan), 0.2004225755740747),
...     (hxacan.MillerIndices(-1, 1, -1, hxacan), 0.19703805522607432),
...     (hxacan.MillerIndices(1, -1, -1, hxacan), 0.19683866569762418),
...     (hxacan.MillerIndices(2, 1, -1, hxacan), 0.20042073834713292),
...     (hxacan.MillerIndices(-2, 2, 0, hxacan), 0.21988965080867295),
...     (hxacan.MillerIndices(-2, -1, -1, hxacan), 0.2007281995185276),
...     (hxacan.MillerIndices(-2, 1, -1, hxacan), 0.20072647044205774),
...     (hxacan.MillerIndices(2, -1, -1, hxacan), 0.20041901523173083),
...     (hxacan.MillerIndices(-2, 1, 1, hxacan), 0.20072463001602556),
...     (hxacan.MillerIndices(0, 2, 0, hxacan), 0.23260015100470804),
...     (hxacan.MillerIndices(1, 1, -1, hxacan), 0.1968395609520126),
...     (hxacan.MillerIndices(-2, 0, 1, hxacan), 0.19266228580032901),
...     (hxacan.MillerIndices(2, 0, 1, hxacan), 0.19235971851036732),
...     (hxacan.MillerIndices(-2, 0, 0, hxacan), 0.1956121055387738),
...     (hxacan.MillerIndices(1, 2, 0, hxacan), 0.2273675403933891),
...     (hxacan.MillerIndices(2, 0, 0, hxacan), 0.1952180443198133),
...     (hxacan.MillerIndices(2, 0, -1, hxacan), 0.19235690444211515),
...     (hxacan.MillerIndices(-2, 0, -1, hxacan), 0.1926651034625492),
...     (hxacan.MillerIndices(-2, -2, 0, hxacan), 0.21988995158315175),
...     (hxacan.MillerIndices(0, -2, 0, hxacan), 0.23260015804584852),
...     (hxacan.MillerIndices(2, 2, 0, hxacan), 0.21953027052221527),
...     (hxacan.MillerIndices(1, 1, 1, hxacan), 0.19683926430195395),
...     (hxacan.MillerIndices(2, -1, 1, hxacan), 0.20042085241706384),
...     (hxacan.MillerIndices(-1, 1, 1, hxacan), 0.19703835382502216),
...     (hxacan.MillerIndices(-1, -1, 1, hxacan), 0.19703925375799888),
...     (hxacan.MillerIndices(-1, 2, 0, hxacan), 0.2276236144133509),
...     (hxacan.MillerIndices(1, -2, 0, hxacan), 0.22736517051933225),
...     (hxacan.MillerIndices(-1, -2, 0, hxacan), 0.2276260079936761),
...     (hxacan.MillerIndices(-2, -1, 1, hxacan), 0.20072635905074157),
...     (hxacan.MillerIndices(1, -1, 1, hxacan), 0.19683836895638027),
...     (hxacan.MillerIndices(2, -2, 0, hxacan), 0.21952998690923803),
...     (hxacan.MillerIndices(-1, -1, -1, hxacan), 0.1970389552502907)
... ]
>>> original_morphology = BFDHMorphology(hxacan)
>>> calculated_morphology = BFDHMorphology.from_growth_rates(hxacan, growth_rates)
>>> print(original_morphology.bounding_box)
(Coordinates(x=-17.548, y=-12.069, z=-22.105), Coordinates(x=17.548, y=12.069, z=22.105))
>>> print(calculated_morphology.bounding_box)
(Coordinates(x=-15.709, y=-18.679, z=-19.865), Coordinates(x=15.677, y=18.679, z=19.865))

You can see that the two morphologies differ markedly.

Symmetry Allowed Morphology

Note

The ccdc.morphology.SymmetryAllowedMorphology features are available only to CSD-Particle users.

BFDH morphology assumes that a facet’s growth rate is inversely proportional to its interplanar spacing. The ccdc.morphology.SymmetryAllowedMorphology class takes a step back from this assumption: it enumerates every facet that is permitted by the crystal’s point-group symmetry within a chosen range of Miller indices, without applying any growth law to rank them. Every enumerated facet is given the same, default relative growth rate, and the resulting habit appears a a multi-faceted polyhedron with more or less spherical shape, depending on the crystal structure that is being analysed.

This morphology to the Symmetry Allowed option of the Growth Explorer’s Morphology tab in Mercury, which is located in the CSD-Particle menu.

As before, we start by reading a crystal structure, this time using the CSD refcode ADIPAC (adipic acid):

>>> from ccdc.io import EntryReader
>>> from ccdc.morphology import SymmetryAllowedMorphology
>>> reader = EntryReader('csd')
>>> crystal = reader.crystal('ADIPAC')

The calculation is controlled by a ccdc.morphology.SymmetryAllowedMorphology.Settings instance:

>>> settings = SymmetryAllowedMorphology.Settings()
>>> print(settings.maximum_hkl_integer)
3
>>> print(settings.minimum_d_spacing)
0.5

• maximum_hkl_integer limits the largest integer used for h, k or l. Increasing it will allow higher order reflections to be considered, and so will generally increase the number of facets found.

• minimum_d_spacing describes the smallest interplanar spacing, in Angstroms, that should be considered.

Both attributes reject non-physical values:

>>> settings.maximum_hkl_integer = 0
Traceback (most recent call last):
    ...
ValueError: maximum_hkl_integer must be a positive integer
>>> settings.minimum_d_spacing = -1.0
Traceback (most recent call last):
    ...
ValueError: minimum_d_spacing must be a positive number

Once we are happy with our settings we can create the morphology:

>>> settings = SymmetryAllowedMorphology.Settings()
>>> morphology = SymmetryAllowedMorphology(crystal, settings)
>>> print(len(morphology.facets)) 
290

Increasing maximum_hkl_integer will typically increase the number of facets returned, since higher order reflections become admissible:

>>> settings.maximum_hkl_integer = 4
>>> wider_morphology = SymmetryAllowedMorphology(crystal, settings)
>>> print(len(wider_morphology.facets)) 
482

Note

It is possible to set the minimum_d_spacing setting, but for the moment this setting has no effect on the calculation and it is reserved for future use. At present, maximum_hkl_integer is the only setting that controls the size of the generated facet list.

Like BFDHMorphology, SymmetryAllowedMorphology is a MorphologyBase, so every property and visualisation technique shown above, such as volume, bounding_box and the wireframe/skin plotting functions, can be reused without any changes:

>>> print(round(morphology.volume, 3)) 
5432.101

A symmetry allowed morphology is generally most useful as a source of candidate Miller indices for further calculations.

Seed Morphology

Note

The ccdc.morphology.SeedMorphology features are available only to CSD-Particle users.

While ccdc.morphology.SymmetryAllowedMorphology and ccdc.morphology.BFDHMorphology generate a morphology from symmetry and geometry alone, ccdc.morphology.SeedMorphology uses the network of intermolecular interactions (Periodic Bond Chains or PBCs) present in the crystal structure to identify the most likely facets following Hartman-Perdok theory. Literature references can be found in Calculating Periodic Bond Chains.

Following Hartman and Perdok’s definition, a crystal plane (or face) can be defined as:

• Kinked (K) if it contains no PBC vector.

• Stepped (S) if it contains one coplanar PBC vector.

• Flat (F) if it contains two or more coplanar PBC vectors.

After calculating PBCs, a morphology is constructed by retaining only the Flat faces, and assigning each of them the same growth rate.

This morphology is the same that can be calculated selecting the Interaction Network option from Growth Explorer in Mercury.

As well as a crystal morphology, SeedMorphology gives access to the full hierarchy of surface sites identified from the PBCs, which allows to explore crystal growth in more detail.

This nested hierarchy is composed of:

• Growth face (GrowthFace) - a crystal face, identified by its Miller indices.

• Layer (Layer) - a layer of molecules parallel to the Miller plane.

• Step edge (StepEdge) - an incomplete edge of a layer, formed by one or more rows of molecules.

• Edge row (EdgeRow) - a single molecular row within a step edge.

• Kink site (KinkSite) - an attachment point where a step edge is completed one molecule at a time.

To summarise, a growth face is made up of layers, a layer of step edges, a step edge of edge rows, and an edge row of kink sites.

Settings

ccdc.morphology.SeedMorphology.Settings is a subclass of ccdc.crystal.Crystal.PeriodicBondChain.Settings, so every setting described in Calculating Periodic Bond Chains is available here too, and controls how the underlying PBCs are calculated. For example:

>>> from ccdc.morphology import SeedMorphology
>>> seed_settings = SeedMorphology.Settings()
>>> seed_settings.force_field = 'UNI'
>>> seed_settings.interaction_energy_threshold = -10.0

Calculating the morphology

We construct the morphology by passing a crystal and, optionally, an instance of the Settings class:

>>> reader = EntryReader('csd')
>>> crystal = reader.crystal('ADIPAC')
>>> seed_morphology = SeedMorphology(crystal)
Hydrogens atom positions for the crystal have been normalised to calculate periodic bond chains.

Note

If no PBCs can be found for the crystal, a UserWarning will be raised and the morphology’s properties, other than identifier, will be unavailable.

During the calculation, several combinations of PBCs are tested, and only the faces that satisfy selection criteria in line with the Hartman and Perdok theory are retained in the final morphology.

These can be accessed through the facets attribute, which returns a tuple of facets which bound the morphology polyhedron. Each of these is an instance of the SeedFacet class, a child class of Facet extended to include the underlying Growth Face.

>>> print(len(seed_morphology.facets)) 
14
>>> facet = seed_morphology.get_facet((0, 0, 2))
>>> print(facet.miller_indices.hkl) 
(0, 0, 2)

Information about all of the tested faces can still be accessed, even if they are not all exposed in the final morphology, through the growth_faces property. Alternatively, we can look for a specific Growth Face using the get_growth_face():

>>> print(len(seed_morphology.growth_faces)) 
22
>>> growth_face = seed_morphology.get_growth_face((0, 0, 2))

If any of the growth_faces was not included in the morphology, it is possible that its network of PBCs did not satisfy our selection criteria. We can investigate the reason for any calculation failure using its status property (see below).

Growth status

Every level of the hierarchy, from a Growth Face to an individual Kink Site, will have an associated status attribute, an instance of the GrowthStatus class.

For example, for a Growth Face:

>>> status = growth_face.status
>>> print(status.is_good) 
True
>>> print(status.status) 
Good

If a feature could not be fully resolved, its status will be Bad or Unfinished, and its reason attribute will describe why:

>>> if status.is_bad:
...     print(status.reason) 

Growth energies

Each of the items in the hierarchy has some associated energies, which are calculated as the sum of the intermolecular interactions that were broken to create each site. When their sign is changed, these energies represent approximate interfacial energies for the site of interest.

The class ccdc.morphology.GrowthEnergies is used to store the energies for each construct, which will be accessible through the relevant energies property.

The energy values that are returned are:

• Terrace energy - includes interactions from the layer towards the ‘outside’ of the crystal, excluding any inter-layer energy.

• Reverse Terrace energy - includes interactions from the step towards the ‘outside’ of the layer, excluding any inter-layer energy.

• Step Edge energy - includes interactions from the step towards the ‘outside’ of the layer, excluding any inter-step energy.

• Reverse Step Edge energy - includes interactions from the step towards the ‘inside’ of the layer, excluding any inter-step energy.

• Kink energy - includes interactions from the step towards the ‘outside’ of the row, excluding any inter-row energy.

• Reverse Kink energy - includes interactions from the kink site towards the ‘inside’ of the row, excluding any inter-row energy.

• Internal energy - includes all interactions within molecules that form part of the site.

• Bulk energy - includes all interactions from the site toward the ‘inside’ of the crystal.

• Interfacial energy - includes all interactions from the site towards the ‘outside’ of the crystal.

A depiction of some of these energies can be found in Figure 2 of

Modeling Step Velocities and Edge Surface Structures during Growth of Non-Centrosymmetric Crystals.

For example, for the (002) layer of ADIPAC:

>>> growth_face = seed_morphology.get_growth_face((0, 0, 2))
>>> layer = growth_face.layers[0]
>>> energies = layer.energies
>>> print(f"""Terrace energy         : {round(energies.terrace, 3)}
... Reverse Terrace energy   : {round(energies.reverse_terrace, 3)}
... Step Edge energy         : {round(energies.step_edge, 1)}
... Reverse Step Edge energy : {round(energies.reverse_step_edge, 1)}
... Internal energy          : {round(energies.internal, 3)}
... Bulk energy              : {round(energies.bulk, 3)}
... Interfacial energy       : {round(energies.interfacial, 3)}""") 
Terrace energy         : -30.874
Reverse Terrace energy : -30.874
Step Edge energy       : 0.0
Reverse Step Edge energy : 0.0
Internal energy        : -56.976
Bulk energy            : -30.874
Interfacial energy     : -30.874

All the energies above have units of kJ/mol.

Not every field of GrowthEnergies is meaningful at every level of the hierarchy, and unused fields are always zero:

Level	Fields that are always zero
Layer	step_edge, reverse_step_edge, kink, reverse_kink
Step edge	kink, reverse_kink
Edge row	kink, reverse_kink
Kink site	internal

Exploring the growth hierarchy

We already have a GrowthFace for the (0, 0, 2) face of ADIPAC. Its properties describe the face as a whole:

>>> print(growth_face.miller_indices) 
(0, 0, 2)
>>> print(round(growth_face.interplanar_spacing, 3)) 
3.446470441341587
>>> print(round(growth_face.average_height_per_layer, 3)) 
3.446470441341587
>>> print(growth_face.num_symmetry_layers) 
2
>>> print(tuple(round(x, 3) for x in growth_face.face_orth_normal)) 
(0.0, 0.0, 1.0)

A face is made up of one or more layers, each a complete, repeating unit of growth perpendicular to the face:

>>> layers = growth_face.layers
>>> print(len(layers)) 
2
>>> layer = layers[0]
>>> print(len(layer.molecules)) 
1
>>> print(layer.step_edge_directions) 
[(0, 1, 0), (0, -1, 0), (1, 0, 0), (-1, 0, 0)]

An incomplete layer grows outwards along a step edge, in the direction perpendicular to one of its constituent PBCs. Cnventionally, a step edge is described by the direction along which it is running, i.e., the direction of its constituent PBCs. This is the step_direction. We can also access the step_growth_direction, which is the direction perpendicular to the step.

>>> step_edge = layer.step_edge_by_direction((0, -1, 0))
>>> print(step_edge.step_direction) 
(0, -1, 0)
>>> print(step_edge.step_growth_direction) 
(1, 0, 0)
>>> print(round(step_edge.average_propagation_length, 3)) 
10.01
>>> edge_rows = step_edge.edge_rows
>>> print(len(edge_rows)) 
1

A step edge is itself made up of edge rows: the symmetry-unique, repeating molecular rows that run along the step edge direction:

>>> edge_row = edge_rows[0]
>>> print(edge_row.step_direction) 
(0, -1, 0)
>>> print(round(edge_row.average_width_per_growth_unit, 3)) 
5.15

Finally, an edge row is completed one molecule at a time at its kink sites: the primary attachment points for new growth units during crystal growth. Each kink site propagates the row in an East or West direction (keeping the direction perpendicular to the step edge as North):

>>> kink_sites = edge_row.kink_sites
>>> print(len(kink_sites)) 
2
>>> kink_site = kink_sites[0]
>>> print(kink_site.kink_direction) 
East
>>> print(kink_site.kink_direction_vector) 
(0, -1, 0)
>>> preceding = kink_site.preceding_molecule
>>> print(preceding.identifier if preceding else None) 
ADIPAC

As with layers and step edges, a kink site’s energies describe the energetics of adding a molecule at that site, which can be used, for example, to estimate its relative rate of growth.

Accessing the molecules of a growth feature

Most of the growth hierarchy describes geometry and energetics, but a few properties and methods return instances of the ccdc.molecule.Molecule class, which can be inspected, or visualised, like any other molecule in the API.

We already saw that a layer’s molecules property tells us how many symmetry-unique molecules make up that layer; we can just as easily inspect the molecules themselves:

>>> layer_molecules = layer.molecules
>>> print([mol.identifier for mol in layer_molecules]) 
['ADIPAC', 'ADIPAC']

Given one of these molecules, get_layer_containing_molecule on a growth face will look up which layer it belongs to, if any:

>>> containing_layer = growth_face.get_layer_containing_molecule(layer_molecules[0])
>>> print(containing_layer is not None) 
True

At the bottom of the hierarchy, a kink site’s preceding_molecule, shown above, returns the molecule immediately behind the kink vacancy in the growth direction: the molecule that a newly attached growth unit would pack against. It returns None if there is no such molecule, for example at the very start of a row.

No other property of the growth hierarchy returns a molecule directly; the remaining geometric and energetic properties describe a growth feature as a whole, rather than any individual molecule within it.