Calculating Periodic Bond Chains

Note

The Periodic Bond Chain features are available only to CSD-Materials and CSD-Enterprise users.

The class ccdc.crystal.Crystal.PeriodicBondChain allows calculation of Periodic Bond Chains (PBC) for a crystal structure.

PBCs were first introduced by Hartman and Perdok as a way to find a connection between crystal structure and crystal morphology, as described in their series of articles published in 1955:

• On the relations between structure and morphology of crystals. I.

• On the relations between structure and morphology of crystals. II.

• On the relations between structure and morphology of crystals. III.

According to Hartman and Perdok’s theory, PBCs can help to identify directions of preferential growth of a crystal by considering uninterrupted chains of strong bonds (or intermolecular interactions) between the crystal’s building units, which can be atoms or, like in our case, molecules.

Each PBC is identified with a vector indicating its direction, as well as information about the molecules that form part of the PBC.

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.

Here, we will see a few examples of how PBCs can be calculated for a crystal structure, and how they can be used to interrogate patterns of intermolecular interactions that are responsible for crystal growth.

As usual, we start by selecting a CSD refcode to initialise our crystal structure:

>>> from ccdc import io
>>> from ccdc.crystal import Crystal
>>> csd_refcode = 'ACSALA01'
>>> csd_reader = io.EntryReader("CSD")
>>> crystal = csd_reader.entry(csd_refcode).crystal

We can then calculate the PBCs for our structure:

>>> periodic_bond_chains = crystal.periodic_bond_chains()
Hydrogens atom positions for the crystal have been normalised to calculate periodic bond chains.
>>> len(periodic_bond_chains)
263

The above calculation used default settings, but we can modify the list of returned PBCs by selecting different settings, as explained in the next section.

We can look at one of the calculated PBCs in more detail by accessing its various attributes:

>>> first_pbc = periodic_bond_chains[0]
>>> print(f"This PBC runs along the {first_pbc.direction} direction and it is formed by {first_pbc.length} interactions")
This PBC runs along the (0, 1, 0) direction and it is formed by 2 interactions
>>> print(f"The strongest intermolecular interaction forming this PBC has a strength of {round(first_pbc.max_strength)} kJ/mol") 
The strongest intermolecular interaction forming this PBC has a strength of -49.5 kJ/mol
>>> print(f"The weakest intermolecular interaction forming this PBC has a strength of {round(first_pbc.min_strength)} kJ/mol") 
The weakest intermolecular interaction forming this PBC has a strength of -23.7 kJ/mol
>>> print(f"The average energy for intermolecular interactions forming this PBC is {round(first_pbc.average_strength)} kJ/mol") 
The average strength of this PBC is -21.5 kJ/mol

The “width” attribute indicates the maximum lateral span (in Angstroms) of the chain. A smaller value will indicate a straighter PBC, whereas a larger value indicates a jagged chain.

If the only interaction forming a chain is that connecting two translation-related molecules, then the chain width will be zero.

>>> other_pbc = periodic_bond_chains[28]
>>> print(f"This PBC is formed by {other.length} interactions and has a width of {round(other_pbc.width, 1)} Angstroms") 
The average strength of this PBC is -21.5 kJ/mol
This PBC is formed by 1 interactions and has a width of 0.0 Angstroms

The last attribute for each PBC is the underlying Synthon Data, i.e. the list of pairwise intermolecular interactions that are forming the chain.

Selecting settings for the calculation of PBCs

The calculation of Periodic Bond Chains for a crystal structure can be controlled using the ccdc.crystal.Crystal.PeriodicBondChain.Settings class.

We first need to instantiate the class

>>> settings = Crystal.PeriodicBondChain.Settings()

Then, after we are happy with our settings, we can calculate the PBCs by passing the settings instance as an argument:

>>> periodic_bond_chains = crystal.periodic_bond_chains(settings=settings) 

We can change the force field used for the calculation (by default, this is set to “CLP”):

>>> settings.force_field = "UNI"

We can change the minimum energy that an interaction needs to have to be included in a PBC

>>> settings.interaction_energy_threshold = -10
>>> periodic_bond_chains = crystal.periodic_bond_chains(settings=settings)
Hydrogens atom positions for the crystal have been normalised to calculate periodic bond chains.
>>> print(f"{len(periodic_bond_chains)} PBCs found using the {settings.force_field} force field with energy threshold {round(settings.interaction_energy_threshold, 1)} kJ/mol")
78 PBCs found using the UNI force field with energy threshold -10 kJ/mol

A larger energy threshold will likely result in fewer PBCs, but it might be desirable if we want to focus only on the strongest interactions in our crystal.

Similarly, we can also reduce the number of calculated PBCs by imposing some constraints on the maximum allowed width (the lateral span of a PBC). To show the difference, we first re-initiate the settings to the default values.

>>> settings = Crystal.PeriodicBondChain.Settings()
>>> settings.chain_width_limit = 4.0
>>> periodic_bond_chains = crystal.periodic_bond_chains(settings=settings)
Hydrogens atom positions for the crystal have been normalised to calculate periodic bond chains.
>>> print(f"{len(periodic_bond_chains)} PBCs found using the {settings.force_field} force field with energy threshold {round(settings.interaction_energy_threshold, 1)} kJ/mol and a chain width limit of {settings.chain_width_limit} Angstroms")
54 PBCs found using the CLP force field with energy threshold -4.2 kJ/mol and a chain width limit of 4.0 Angstroms

Or we can change the “minimum forward progress limit”. This setting controls the minimum step distance in the direction of a PBC that an interaction needs to have in order to be included. A negative value allows for interactions in the chain to connect molecules in a direction opposite to the overall PBC direction, whereas a positive value will prevent this. To show the difference, we first re-initiate the settings to the default values.

Note that this setting also needs to be enabled before the calculation.

>>> settings = Crystal.PeriodicBondChain.Settings()
>>> settings.min_forward_progress_limit = 4.0
>>> settings.enable_min_forward_progress_limit = True
>>> periodic_bond_chains = crystal.periodic_bond_chains(settings=settings) 
>>> print(f"{len(periodic_bond_chains)} PBCs found using the {settings.force_field} force field with energy threshold {round(settings.interaction_energy_threshold, 1)} kJ/mol and a minimum forward progress limit of {settings.min_forward_progress_limit} Angstroms") 
115 PBCs found using the CLP force field with energy threshold -4.2 kJ/mol and a minimum forward progress limit of 4.0 Angstroms

The remaining settings are:

• Crystal.PeriodicBondChain.Settings.interaction_distance_limit

  controls the distance limit (in Angstroms) for interactions that will be included in the PBC calculation (this is 5.0 Angstroms by default). A larger value will find more PBCs.

• Crystal.PeriodicBondChain.Settings.normalise_hydrogens

  controls wheter the positions of hydrogen atoms are normalised using neutron diffraction data before calculating PBCs (this is True by default)

The PBC Synthon Data class

Each intermolecular interaction in the pbc is represented as an instance of the ccdc.crystal.Crystal.PeriodicBondChain.PBCSynthonData class.

This class has attributes that describe the intermolecular interaction and the molecules that are involved in the latter. A “from” molecule is the molecule from which the interaction generates, while a “to” molecule is the molecule receiving the interaction. The assignment of the “from” and “to” properties is based on some conventions, and we can easily invert the assignment without affecting the PBC, as we will see later.

>>> first_synthon = periodic_bond_chains[0].synthon_data[0]
>>> first_synthon
SynthonData(0, 'x,y,z', (0,0,0), 0, '-x,-y,-z', (-1,1,-1), -49.541, 7.777)
>>> print(f"This interaction has strength {round(first_synthon.strength, 1)} kJ/mol and distance {round(first_synthon.distance, 3)} Angstroms")
This interaction has strength -49.5 kJ/mol and distance 7.777 Angstroms

Each of the molecules involved in the interaction is identified by:

• unit cell translation: a tuple((int, int, int)) representing the unit cell containing the molecule.

  (0, 0, 0) indicates the original unit cell, while (1, 1, 1), for example, indicates the unit cell translated by one a axis vector, one b axis vector and one c axis vector.

• symmetry operation: the symmetry operation that relates the molecule of interest and the original molecule in the crystal’s asymmetric unit

We can access these attributes like this:

>>> first_synthon.from_mol_translation
(0, 0, 0)
>>> first_synthon.from_mol_symmop 
x,y,z
>>> first_synthon.to_mol_translation
(-1, 1, -1)
>>> first_synthon.to_mol_symmop 
-x,-y,-z

So, for this synthon, we see that the interaction with strength of -49.5 kJ/mol and distance 7.777 Angstroms is connecting a molecule in the asymmetric unit (unit cell (0,0,0) and symmop “x,y,z”) with an inversion-related molecule (symmop is “-x,-y-z”) in the cell translated by -1a, 1b and -1c.

Performing Operations on PBCs

As mentioned previously, each of the calculated PBCs follows some conventions to determine the starting molecule and the chain direction. We can look at other symmetry-related PBCs in our crystal by performing some operations on the original PBCs.

The first operation is inversion.

This operation inverts the direction of the whole chain, but without changing the interactions in the chain. Simply, the last interaction in the original chain will become the last in the inverted case. The new inverted PBC will be a copy of the original one.

>>> original_pbc = periodic_bond_chains[0]
>>> inverted_pbc = original_pbc.reverse()
>>> print(f"Direction of original PBC: {original_pbc.direction}. Direction of inverted PBC: {inverted_pbc.direction}")
Direction of original PBC: (0, 1, 0). Direction of inverted PBC: (0, -1, 0)
>>> original_pbc.synthon_data[0]
SynthonData(0, 'x,y,z', (0,0,0), 0, '-x,-y,-z', (-1,1,-1), -49.541, 7.777)
>>> inverted_pbc.synthon_data[0]
SynthonData(0, 'x,y,z', (0,1,0), 0, '-x,-y,-z', (-1,1,-1), -23.669, 5.568)

The next possible operation is translation.

This operation will return a copy of the PBC translated to another unit cell, specified by a translation component. In the example below, (1,1,1) represent a translation of one a axis vector, one b axis vector and one c axis vector.

>>> original_pbc = periodic_bond_chains[0]
>>> translated_pbc = original_pbc.translate((1, 1, 1))
>>> print(f"Direction of original PBC: {original_pbc.direction}. Direction of translated PBC: {translated_pbc.direction}")
Direction of original PBC: (0, 1, 0). Direction of translated PBC: (0, 1, 0)
>>> original_pbc.synthon_data[0]
SynthonData(0, 'x,y,z', (0,0,0), 0, '-x,-y,-z', (-1,1,-1), -49.541, 7.777)
>>> translated_pbc.synthon_data[0]
SynthonData(0, 'x,y,z', (1,1,1), 0, '-x,-y,-z', (0,2,0), -49.541, 7.777)

The final operation is rotation.

This is the least intuitive of the operations available, but having access to the synthon data helps us to keep track of things. The reason why this operation might not be intuitive is because we are not performing a rotation in the conventional sense (i.e. around an axis), but we are rotating over the molecules forming a PBC.

So, if a PBC made of molecules A and B can be represented as:

A(1)-B(1)-A(2)

a rotation of one step will return the equivalent PBC starting on B(1):

B(1)-A(2)-B(2)

This rotation is specified using an integer number of steps, representing the number of interactions that are “skipped”. We can also use a negative value, which will result in navigating the PBC in the opposite direction.

As an example using a single step for the rotation:

>>> original_pbc = periodic_bond_chains[0]
>>> rotated_pbc = original_pbc.rotate(steps=1)
>>> print(f"Direction of original PBC: {original_pbc.direction}. Direction of rotated PBC: {rotated_pbc.direction}")
Direction of original PBC: (0, 1, 0). Direction of rotated PBC: (0, 1, 0)
>>> print("Synthons of the original PBC")
Synthons of the original PBC
>>> original_pbc.synthon_data[0]
SynthonData(0, 'x,y,z', (0,0,0), 0, '-x,-y,-z', (-1,1,-1), -49.541, 7.777)
>>> original_pbc.synthon_data[1]
SynthonData(0, '-x,-y,-z', (-1,1,-1), 0, 'x,y,z', (0,1,0), -23.669, 5.568)
>>> print("Synthons of the rotated PBC")
Synthons of the rotated PBC
>>> rotated_pbc.synthon_data[0]
SynthonData(0, '-x,-y,-z', (-1,1,-1), 0, 'x,y,z', (0,1,0), -23.669, 5.568)
>>> rotated_pbc.synthon_data[1]
SynthonData(0, 'x,y,z', (0,1,0), 0, '-x,-y,-z', (-1,2,-1), -49.541, 7.777)

Let’s see another example using a negative step this time. Here, we will start two steps backwards, and since the PBC is formed of two interactions, we will end up on the original asymmetric unit molecule (with symmop “x,y,z” and translation (0,0,0)).

>>> original_pbc = periodic_bond_chains[0]
>>> rotated_pbc = original_pbc.rotate(steps=-2)
>>> print(f"Direction of original PBC: {original_pbc.direction}. Direction of rotated PBC: {rotated_pbc.direction}")
Direction of original PBC: (0, 1, 0). Direction of rotated PBC: (0, 1, 0)
>>> print("Synthons of the original PBC")
Synthons of the original PBC
>>> original_pbc.synthon_data[0]
SynthonData(0, 'x,y,z', (0,0,0), 0, '-x,-y,-z', (-1,1,-1), -49.541, 7.777)
>>> original_pbc.synthon_data[1]
SynthonData(0, '-x,-y,-z', (-1,1,-1), 0, 'x,y,z', (0,1,0), -23.669, 5.568)
>>> print("Synthons of the rotated PBC")
Synthons of the rotated PBC
>>> rotated_pbc.synthon_data[0]
SynthonData(0, 'x,y,z', (0,-1,0), 0, '-x,-y,-z', (-1,0,-1), -49.541, 7.777)
>>> rotated_pbc.synthon_data[1]
SynthonData(0, '-x,-y,-z', (-1,0,-1), 0, 'x,y,z', (0,0,0), -23.669, 5.568)

Note that the overall direction for the PBC remains unchanged.