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

  1. Filippini and A. Gavezzotti, Acta Cryst. B49 (1993) 868-880

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

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]

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.