VisualHabit¶
The class ccdc.morphology.VisualHabit provides a simple programmatic interface to the software VisualHabit.
VisualHabit uses the atom-atom pair-wise summation to determine the intermolecular interactions in a molecular crystal.
Summing all the atom-atom interactions across all molecules up to a limiting radius yields the lattice energy.
Calculating the interactions along specific crystallographic directions allows the slice and attachment energies to be calculated.
The attachment energy is a measure of the relative growth rate along specific directions and consequently the growth morphology can be modelled. Please see VisualHabit95 — a program for predicting the morphology of molecular crystals as a function of the growth environment for details of the VisualHabit software.
In order to make use of VisualHabit, we first instantiate the VisualHabit class:
>>> from ccdc.morphology import VisualHabit
>>> visualhabit_settings = VisualHabit.Settings()
he VisualHabit class, as is usual with a complex class in the CSD Python API, has a nested :class:ccdc.morphology.VisualHabit.Settings class. Here it is used to specify a forcefield to use - by default the Dreiding II forcefield (S.L. Mayo, B.D. Olafson and W.A. Goddard III, J. Phys. Chem. 1990 (94) 8897) which contains parameters for the atom types to be found in organic molecules as well as a small set of metal atoms. There are other more specialised forcefields available, such as: Momany (F.A. Momany, L.M. Carruthers, R.F. McGuire and H.A. Scheraga, J. Phys. Chem. 78 (1974) 1595) for a purely organic set of atom types; Gavezotti (G. Filippini and A. Gavezzotti, Acta Cryst. B49 (1993) 868-880) again for organic crystals.
The potential can be set in the Settings class by assigning any unambiguous prefix of the name of the forcefield to the ccdc.morphology.VisualHabit.Settings.potential attribute.
>>> visualhabit_settings.potential = 'dreidingII'
would set the potential to ‘DreidingII’.
The full set of available potentials is:
Name |
Reference |
|---|---|
dreidingII |
S.L. Mayo, B.D. Olafson and W.A. Goddard III, J. Phys. Chem. 1990 (94) 8897 |
gavezzotti |
|
momany |
F.A. Momany, L.M. Carruthers, R.F. McGuire and H.A. Scheraga, J. Phys. Chem. 78 (1974) 1595 |
If the forcefield selected is inapplicable to the crystal for which we wish to calculate lattice energies, a RuntimeError will be raised.
The other parameter of the VisualHabit calculation is the convergence limit. This is a distance in Angstroms from the centre of the crystal packing shell to calculate intermolecular energies. The default value is 30.0, which is reasonable in most cases, but may be set to any desired value.
Once the settings have been determined, the calculation may be performed on any crystal.
In this example we will calculate the lattice energy of the drug Nabumetone (XOCXUI02)
First we load the crystal structure XOCXUI02 from the CSD in the usual way:
>>> from ccdc import io
>>> csd = io.EntryReader('CSD')
>>> nabumetone = csd.crystal('XOCXUI02')
Next we calculate the lattice energy using VisualHabit:
>>> results = VisualHabit(settings=visualhabit_settings).calculate(nabumetone)
This will return an instance of ccdc.morphology.VisualHabit.Results from which various energy terms, and a predicted morphology may be extracted.
First we may access the Lattice energy from the results by examining results.lattice_energy which provides an access to the total lattice energy and its components. These energies are those of the whole crystal’s lattice:
>>> lattice_energy = results.lattice_energy
>>> print(f"""Total Latt Energy : {round(lattice_energy.total, 3)} kJ/mol
... Electrostatic Energy : {round(lattice_energy.electrostatic, 3)} kJ/mol
... H-Bond Energy : {round(lattice_energy.h_bond, 3)} kJ/mol
... H-Bond Attraction Energy : {round(lattice_energy.h_bond_attraction, 3)} kJ/mol
... H-Bond Repulsion Energy : {round(lattice_energy.h_bond_repulsion, 3)} kJ/mol
... van deer Waals Energy : {round(lattice_energy.vdw, 3)} kJ/mol
... van deer Waals Attraction Energy : {round(lattice_energy.vdw_attraction, 3)} kJ/mol
... van deer Waals Repulsion Energy : {round(lattice_energy.vdw_repulsion, 3)} kJ/mol""")
Total Latt Energy : -123.084 kJ/mol
Electrostatic Energy : -11.22 kJ/mol
H-Bond Energy : 0.0 kJ/mol
H-Bond Attraction Energy : 0.0 kJ/mol
H-Bond Repulsion Energy : 0.0 kJ/mol
van deer Waals Energy : -111.864 kJ/mol
van deer Waals Attraction Energy : -218.141 kJ/mol
van deer Waals Repulsion Energy : 106.278 kJ/mol
We can access the list of synthons (the intermolecular interactions) that make up the lattice energy using the synthons property.
>>> synthons = results.synthons
>>> print(f"Number of synthons {len(synthons)}")
Number of synthons 365
Each synthon contains the ccdc.molecule.Molecule pairs involved in the itermolecular interactions, the distances between their centroids (center of geometry) and computed engergies.
As an example we can look at the first two synthons in the list and check how many acceptors exist in each molecule. Of course we expect the number of acceptors to be the same in each molecule, as this is a single-component system.
>>> for syn in synthons[:2]:
... print(f"Synthon distance : {round(syn.distance,3)}")
... print(f"Total Energy: {round(syn.total,3)} kJ/mol.")
... print(f"HB Energy: {syn.h_bond} kJ/mol.")
... print(f"Electrostatic Energy: {round(syn.electrostatic,3)} kJ/mol.")
... print(f"Van der Waals Energy: {round(syn.vdw,3)} kJ/mol.")
... print(f"Number of Acceptors in Mol_1 : {len([atom for atom in syn.molecule1.atoms if atom.is_acceptor])}")
... print(f"Number of Acceptors in Mol_2 : {len([atom for atom in syn.molecule2.atoms if atom.is_acceptor])}")
Synthon distance : 4.652
Total Energy: -19.718 kJ/mol.
HB Energy: 0.0 kJ/mol.
Electrostatic Energy: -3.521 kJ/mol.
Van der Waals Energy: -16.197 kJ/mol.
Number of Acceptors in Mol_1 : 2
Number of Acceptors in Mol_2 : 2
Synthon distance : 4.652
Total Energy: -19.718 kJ/mol.
HB Energy: 0.0 kJ/mol.
Electrostatic Energy: -3.521 kJ/mol.
Van der Waals Energy: -16.197 kJ/mol.
Number of Acceptors in Mol_1 : 2
Number of Acceptors in Mol_2 : 2
We can also access the symmetry operator and tranlation vector between Mol_1 to Mol_2:
>>> syn_1 = synthons[0]
>>> print(f"Symmetry operator : {syn_1.symmetry_operator}")
Symmetry operator : x,1/2-y,1/2+z
>>> print(f"Unit Cell translation: {syn_1.translation_vector}")
Unit Cell translation: (0, 0, 0)
From the results class, a ccdc.morphology.VisualHabitMorphology may be obtained.
>>> morph = results.morphology
An interesting scientific use-case for the VisualHabit API module is the ability to assess the calculated strengths of different types of interaction on the faces of a crystal. For example, crystal faces with strong hydrogen bonding interactions may experience strong interactions with water, potentially slowing the growth of that face in this solvent. Faces with strong electrostatic interactions may aggregate, or stick to machinery in an industrial production line. Understanding the types of interactions likely to form at the largest crystal faces may offer possible insight into crystallisation and solid-state behaviour.
To analyze the interactions at each face of a crystal we first need to list the predicted faces of the morphology.
>>> faces = results.morphology.facets
It can be useful to look at the Miller index of each face, and also which faces are dominant (largest) in the morphology. This can be done by sorting a list of faces by their relative area.
First we create a nested list of miller indices and their relative areas:
>>> face_list = [[face.miller_indices.hkl, face.area] for face in faces]
We can also access the offset value (shift along plane normal), which is optimised to achieve the most stable surface termination. This value can be used for further analysis within the ccdc.particle.Surface class.
>>> hkl_and_offsets = [[face.miller_indices.hkl, face.offset] for face in faces]
>>> for indices, offset in hkl_and_offsets:
... print(f"hkl: {indices} | offset: {round(offset,2)}")
hkl: (1, 0, 0) | offset: 0.0
hkl: (-1, 0, 0) | offset: 0.0
hkl: (1, 1, 0) | offset: 1.15
hkl: (1, -1, 0) | offset: 1.15
hkl: (-1, 1, 0) | offset: 1.15
hkl: (-1, -1, 0) | offset: 1.15
hkl: (0, 1, 1) | offset: 0.47
hkl: (0, 1, -1) | offset: 0.47
hkl: (0, -1, 1) | offset: 0.47
hkl: (0, -1, -1) | offset: 0.47
hkl: (1, 1, 1) | offset: 0.46
hkl: (1, -1, 1) | offset: 0.46
hkl: (-1, 1, -1) | offset: 0.46
hkl: (-1, -1, -1) | offset: 0.46
hkl: (1, 0, -2) | offset: -0.78
hkl: (-1, 0, 2) | offset: -0.78
hkl: (0, 0, 2) | offset: 0.0
hkl: (0, 0, -2) | offset: 0.0
hkl: (1, 0, 2) | offset: 0.39
hkl: (-1, 0, -2) | offset: 0.39
Then we sort the list of lists by the area value, using a lambda function to specify the second element (area):
>>> face_list = sorted(face_list, key = lambda x: x[1], reverse=True)
>>> for indices, area in face_list:
... print(f"hkl: {indices} | Area: {round(area, 3)}")
hkl: (-1, 0, 0) | Area: 1850.748
hkl: (1, 0, 0) | Area: 1850.748
hkl: (0, -1, -1) | Area: 151.395
hkl: (0, -1, 1) | Area: 151.395
hkl: (0, 1, -1) | Area: 151.395
hkl: (0, 1, 1) | Area: 151.395
hkl: (1, -1, 0) | Area: 88.699
hkl: (1, 1, 0) | Area: 88.699
hkl: (-1, 1, 0) | Area: 88.699
hkl: (-1, -1, 0) | Area: 88.699
hkl: (0, 0, 2) | Area: 87.237
hkl: (0, 0, -2) | Area: 87.237
hkl: (1, 0, -2) | Area: 58.496
hkl: (-1, 0, 2) | Area: 58.496
hkl: (1, 0, 2) | Area: 3.847
hkl: (-1, 0, -2) | Area: 3.847
hkl: (1, 1, 1) | Area: 3.201
hkl: (-1, 1, -1) | Area: 3.201
hkl: (-1, -1, -1) | Area: 3.201
hkl: (1, -1, 1) | Area: 3.201
We can also analyze the attachment energy of each facet using the method:
>>> face_list = [[face.miller_indices.hkl, face.attachment_energy] for face in faces]
>>> face_list = sorted(face_list, key = lambda x: x[1], reverse=True)
>>> for indices, e_att in face_list:
... print(f"hkl: {indices} | E_Att: {round(e_att, 3)}")
hkl: (1, 0, 0) | E_Att: -12.83
hkl: (-1, 0, 0) | E_Att: -12.83
hkl: (1, 1, 0) | E_Att: -68.106
hkl: (1, -1, 0) | E_Att: -68.106
hkl: (-1, 1, 0) | E_Att: -68.106
hkl: (-1, -1, 0) | E_Att: -68.106
hkl: (0, 1, 1) | E_Att: -75.96
hkl: (0, 1, -1) | E_Att: -75.96
hkl: (0, -1, 1) | E_Att: -75.96
hkl: (0, -1, -1) | E_Att: -75.96
hkl: (1, 1, 1) | E_Att: -76.689
hkl: (1, -1, 1) | E_Att: -76.689
hkl: (-1, 1, -1) | E_Att: -76.689
hkl: (-1, -1, -1) | E_Att: -76.689
hkl: (1, 0, -2) | E_Att: -86.341
hkl: (-1, 0, 2) | E_Att: -86.341
hkl: (0, 0, 2) | E_Att: -86.428
hkl: (0, 0, -2) | E_Att: -86.428
hkl: (1, 0, 2) | E_Att: -86.875
hkl: (-1, 0, -2) | E_Att: -86.875
Comparing the values of the attachment energy and area it becomes apparent that the largest facet exhibited on the morphology is the one with the least negative attachment energy due to the lack of intermolecular interactions present at the surface.
We can also investigate Miller indices that are explicitly defined. These could be orientations that are identified as slip planes from ccdc.particle.SlipPlanes or defined from other calculations.
We provide a hkl, either as a ccdc.crystal.Crystal.MillerIndices object or tuple, we will use the latter in this example.
>>> new_facet = results.facet_energy(miller_indices=(0,1,0))
This will utilise the lattice energy already computed and return the most energetically stable slice. Note - As this facet does not appear on the convex hull the usual area and geometrical descriptors relating to the Convex Hull are not valid. We can now examine the attachment energy and offset
>>> print(f"hkl of the new facet {new_facet.miller_indices.hkl} with an offset of {round(new_facet.offset,3)} A.")
hkl of the new facet (0, 1, 0) with an offset of -0.588 A.
>>> print(f"Attachment energy is {round(new_facet.attachment_energy,3)} kJ/mol.")
Attachment energy is -66.422 kJ/mol.