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.BFDHMorphology.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.