Crystal examples

Generating molecular assemblies

These examples take a crystal and generate a molecular representation based on various criteria.

An assembly based on distance

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2015-06-17: created by the Cambridge Crystallographic Data Centre
#

"""
Write out heaviest component and molecular shell for a set of CSD structures.

This script takes an input a gcd file and an output directory. For each CSD
entry the script then identifies the heaviest component and writes it to the
output directory. The molecules in the molecular shell around the heaviest
component are then identified and written out to a separate file in the output
directory.

For more help use the -h argument.

"""

import argparse
import os
from ccdc import io

##############################################################################

def main(gcd_file, out_dir):
    """
    Run the main program.
    """
    crystal_reader = io.CrystalReader(gcd_file, format='identifiers')
    for crystal in crystal_reader:

        # Get the atoms from the heaviest component in the crystal. These will
        # be used as the selection around which we create the molecular shell.
        heaviest_component = crystal.molecule.heaviest_component
        selection = heaviest_component.atoms

        # Write out the central molecule
        out_file_name = '%s.mol2' % crystal.molecule.identifier
        out_file_path = os.path.join(out_dir, out_file_name)
        with io.MoleculeWriter(out_file_path) as mol_writer:
            mol_writer.write(heaviest_component)

        # Generate the molecular shell around the atoms in the selection.
        mol_shell = crystal.molecular_shell(distance_type="VdW",
                                            distance_range=(0.0, 0.5),
                                            atom_selection=selection)
        # Write out the molecular shell
        out_file_name = '%s_shell.mol2' % crystal.molecule.identifier
        out_file_path = os.path.join(out_dir, out_file_name)
        with io.MoleculeWriter(out_file_path) as mol_writer:
            mol_writer.write(mol_shell)

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('gcd_file', type=str, help="Input GCD file")
    parser.add_argument('out_dir', type=str, help="Output directory")
    args = parser.parse_args()

    if not os.path.isdir(args.out_dir):
        os.mkdir(args.out_dir)
    if not os.path.isfile(args.gcd_file):
        parser.error('%s is not a file.' % args.gcd_file)

    main(args.gcd_file, args.out_dir)

An assembly based on graph sets

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2016-03-23: created by the Cambridge Crystallographic Data Centre
#
'''
    graph_set_closure.py    -   expand a crystal such that all HBonds of the crystal's graph sets
    are included.

    The graph sets of the crystal will be calculated, then the crystal expanded to include
    all of the hydrogen bonds making up the graph sets.
'''
############################################################################

import os
import argparse

from ccdc.descriptors import CrystalDescriptors
from ccdc.io import CrystalReader, MoleculeWriter
from ccdc.molecule import Molecule

############################################################################

class Runner(argparse.ArgumentParser):
    '''Does the whole job.'''
    def __init__(self):
        '''Defines arguments and parses them.'''
        super(self.__class__, self).__init__(
            description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
        )
        self.add_argument(
            'refcode_or_files', nargs='+',
            help='Refcode or crystal file name.'
        )
        self.add_argument(
            '-o', '--output-file', '--output_file',
            help='Output molecule file.  If not given the molecule will be written to stdout in mol2 format.'
        )

        self.args = self.parse_args()

    @property
    def csd(self):
        '''Lazily create a CSD reader if necessary.'''
        if not hasattr(self, '_csd'):
            self._csd = CrystalReader('csd')
        return self._csd

    @property
    def searcher(self):
        '''Lazily create a GraphSetSearch instance.'''
        if not hasattr(self, '_searcher'):
            self._searcher = CrystalDescriptors.GraphSetSearch()
        return self._searcher

    def expand(self, crystal):
        '''Expand the crystal by the symmetry operators of the GraphSet HBonds.'''
        gsets = self.searcher.search(crystal)
        symms = {
            s
            for gs in gsets
            for hb in gs.hbonds
            for s in hb.symmetry_operators
            if s
        }
        mols = [
            crystal.symmetric_molecule(s)
            for s in sorted(symms)
        ]
        mol = Molecule(crystal.identifier)
        for m in mols:
            mol.add_molecule(m)
        return mol

    def run(self):
        '''Runs over all refcodes and files.'''
        if self.args.output_file:
            self.writer = MoleculeWriter(self.args.output_file)
        else:
            self.writer = MoleculeWriter('stdout', format='mol2')
        for r_or_f in self.args.refcode_or_files:
            if os.path.exists(r_or_f):
                reader = CrystalReader(r_or_f)
                for c in reader:
                    self.writer.write(self.expand(c))
            else:
                c = self.csd.crystal(r_or_f)
                self.writer.write(self.expand(c))
        self.writer.close()

############################################################################

if __name__ == '__main__':
    Runner().run()

Other molecular expansions

This example constructs various molecular expansions of a crystal based on packing, slicing, hbond and contact networks, and polymer expansion:

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2017-07-17: created by the Cambridge Crystallographic Data Centre
#

"""
    crystal_expansion.py    -   show the various molecular expansions of a crystal.

    This script creates various multi-molecular structures from a crystal.

    The structures will be written to a multi-mol2 file.
"""

import argparse
from ccdc import io
from ccdc.descriptors import MolecularDescriptors
import os


class Runner(argparse.ArgumentParser):
    """Defines arguments, parses the command line and runs the job."""

    def __init__(self):
        """Defines the arguments."""
        super(self.__class__, self).__init__(description=__doc__,
                                             formatter_class=argparse.ArgumentDefaultsHelpFormatter)
        self.add_argument(
            'refcodes', nargs='+', help='CSD Refcodes to expand.'
        )
        self.add_argument(
            '-o', '--output-directory', '--output_directory', default='.',
            help='Optional output directory'
        )

    def run(self):
        """Runs the task."""
        args = self.parse_args()

        if not os.path.exists(args.output_directory):
            os.makedirs(args.output_directory)

        csd = io.EntryReader('csd')

        for refcode in args.refcodes:
            crystal = csd.crystal(refcode)

            molecular_shell = crystal.molecular_shell(distance_range=(0.0, 7.0))
            molecular_shell.identifier = '%s_molecular_shell' % refcode

            packing_shell = crystal.packing_shell(30)
            packing_shell.identifier = '%s_packing_shell' % refcode

            hbond_network = crystal.hbond_network(repetitions=3)
            hbond_network.identifier = '%s_hbond_network' % refcode

            contact_network = crystal.contact_network(repetitions=2)
            contact_network.identifier = '%s_contact_network' % refcode

            packing = crystal.packing(((-1, -1, -1), (0.5, 0.5, 0.5)))
            packing.identifier = '%s_packing' % refcode

            slicing = crystal.slicing(
                MolecularDescriptors.atom_plane(*[a for a in crystal.molecule.heavy_atoms])
            )
            slicing.identifier = '%s_slicing' % refcode

            polymer_expansion = crystal.polymer_expansion(repetitions=12)
            polymer_expansion.identifier = '%s_polymer' % refcode

            output_file = os.path.join(args.output_directory, '%s.mol2' % refcode)
            with io.MoleculeWriter(output_file) as writer:
                writer.write(molecular_shell)
                writer.write(packing_shell)
                writer.write(hbond_network)
                writer.write(contact_network)
                writer.write(packing)
                writer.write(slicing)
                writer.write(polymer_expansion)


if __name__ == '__main__':
    Runner().run()

Analysis of HBond dimensionality

This example uses the hbond_network expansions and a principle axes aligned molecule to determine in how many dimensions a crystal’s HBonds are expanding.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2017-09-04: created by the Cambridge Crystallographic Data Centre
#
'''
    hbond_dimensionality.py - establish the dimensionality of a crystal's hydrogen bonded network.

    Expands the hbond network of a crystal, and sees which dimenions are growing.  For example:
        ACSALA (aspirin) shows no dimension grows, i.e the network forms a ring
        YIGPIO02 (Ritonavir) grows in one dimension
        HXACAN (Paracetamol) grows in two dimensiont
        DMANTL08 (Mannitol) grows in all three dimensions.
'''
###########################################################################

import os
import argparse

from ccdc import io
from ccdc.descriptors import MolecularDescriptors

csd = io.EntryReader('csd')

###########################################################################

class Runner(argparse.ArgumentParser):
    '''Runs the calculation.'''
    def __init__(self):
        super(self.__class__, self).__init__(description=__doc__)
        self.add_argument('crystals', nargs='+',
            help='Refcodes or crystal files to analyese'
        )
        self.args = self.parse_args()

    def run(self):
        print('Identifier\tNetwork dimensionality')
        for a in self.args.crystals:
            if not os.path.exists(a):
                c = csd.crystal(a)
                self.analyse(c)
            else:
                for c in io.CrystalReader(a):
                    self.analyse(c)

    def analyse(self, crystal):
        '''Calculates the dimensionality of the crystal.'''
        m = crystal.molecule
        siteless = not m.all_atoms_have_sites
        if siteless:
            m.add_hydrogens()
            crystal.molecule = m
        shells = [
            crystal.hbond_network(repetitions=2), crystal.hbond_network(repetitions=5)
        ]
        boxes = [MolecularDescriptors.PrincipleAxesAlignedBox(s) for s in shells]
        vecs = [
            [b.x_vector.length, b.y_vector.length, b.z_vector.length]
            for b in boxes
        ]
        ratios = [vecs[-1][i]/vecs[0][i] for i in range(3)]
        thresh = 1.1
        ndims = sum(r > thresh for r in ratios)
        #print(ratios)
        #print(ndims)
        print('%8s\t%d' % (crystal.identifier, ndims))


###########################################################################

if __name__ == '__main__':
    Runner().run()

Descriptor HTML report

Note that this script makes use of functionality from the cookbook utility module.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
# 2015-06-17: created by the Cambridge Crystallographic Data Centre
#

"""
Generate a HTML table of molecular and crystallographic properties for
structures in the CSD.

This includes, amongst more generic ones, point group analysis descriptors.

For more help use the -h argument.

"""

import sys
import argparse
import os.path
import collections
import codecs

this_dir = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(1, this_dir)
import utilities
sys.path.remove(this_dir)

from ccdc import io
from ccdc.descriptors import MolecularDescriptors

##############################################################################

# The headers for the HTML table
HEADERS = (
    'REFCODE',
    'Identifier',
    'Hydrate',
    'Formula',
    'Space Gp. Symbol',
    'Z Prime',
    'Unique Chem. Units',
    'Mol. Wt.',
    'Donors',
    'Acceptors',
    'Rotatable bonds',
    'Point Gp. Symbol',
    'Point Gp. Order',
    'Point Gp. Description',
    'Atom counts',
)

##############################################################################

def get_num_unique_structures(molecule):
    "Return the number of unique structures in a molecule."
    d = collections.defaultdict(int)
    for m in molecule.components:
        smi = m.smiles
        if smi is None:
            smi = 'Unknown bond in SMILES'
        d[smi] += 1
    return len(d)

def get_atom_counts(molecule):
    "Return a list of strings of atom counts."
    atom_cts = collections.defaultdict(int)
    for atom in molecule.atoms:
        atom_cts[atom.atomic_symbol] += 1
    l = sorted([(v, k) for k, v in atom_cts.items()])
    l.reverse()
    return ['%s: %d' % (k, v) for v, k in l]

##############################################################################

def get_row_data(entry):
    "Return a tuple of molecular and crystallographic properties."

    # Extract the molecule (which may have several components) and the heaviest
    # component from the crystal structure. The heaviest component is the
    # molecule of interest.
    crystal = entry.crystal
    mol = crystal.molecule
    mol_of_interest = mol.heaviest_component

    # Find out how many unique molecules are present in the crystal structure
    num_unique_structures = get_num_unique_structures(mol)

    # Ensure all hydrogen atoms are present, for donor calculations
    mol_of_interest.add_hydrogens()

    # Counts of molecular properties
    donors = len([at for at in mol_of_interest.atoms if at.is_donor])
    accs   = len([at for at in mol_of_interest.atoms if at.is_acceptor])
    rots   = len([b  for b  in mol_of_interest.bonds if b.is_rotatable])

    # Collect a list of atoms present (which could be taken from the formula)
    atom_counts = get_atom_counts(mol_of_interest)

    #   Another datum requested
    hydrate = 'hydrate' in entry.chemical_name.lower()

    # Molecular weight (before stripping hydrogen atoms!)
    mol_wt = mol_of_interest.molecular_weight

    #   Strip hydrogens for point group calculations
    mol_of_interest.remove_hydrogens()

    #   Returns a triple (order, symbol, description)
    pt_group_info = MolecularDescriptors.point_group_analysis(mol_of_interest)

    #   Now we have all the information
    return (
        mol.identifier,
        entry.chemical_name,
        hydrate,
        entry.formula,
        crystal.spacegroup_symbol,
        crystal.z_prime,
        num_unique_structures,
        mol_wt,
        donors,
        accs,
        rots,
        pt_group_info[0],
        pt_group_info[1],
        pt_group_info[2],
        ' '.join(atom_counts)
    )

def main(gcd_file, out_file):
    "Run the main program."

    # Collect the data of interest in a list.
    rows = []
    with io.EntryReader(gcd_file, format='identifiers') as entry_reader:
        for entry in entry_reader:
            rows.append(get_row_data(entry))

    # Write out the report as a HTML file.
    with codecs.open(out_file, 'w', encoding='utf8') as writer:
        writer.write('\n'.join([
            '<HTML>',
            '<HEAD>',
            '<TITLE>Point group analysis</TITLE>',
            '</HEAD>',
            '<BODY>',
            '<H1>Point group analysis</H1>\n'
        ]))
        utilities.write_html_table(HEADERS, rows, stream=writer, border=1)
        writer.write('</BODY></HTML>\n')

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('gcd_file', type=str, help="Input GCD file")
    parser.add_argument('out_file', type=str, help="Output HTML file")
    args = parser.parse_args()

    if not os.path.isfile(args.gcd_file):
        parser.error('%s is not a file.' % args.gcd_file)

    main(args.gcd_file, args.out_file)

Analysing crystal morphologies

This example creates BFDH morphologies of CSD structures and analyses the aspect ratios of them to predict whether they will from needles, plates, laths or blocks. For example, for the polymorphs of the anti-inflammatory XOCXUI01 and XOCXUI02:

Classification of morphologies of XOCXUI01 and XOCXUI02

Classification of XOCXUI01 and XOCXUI02

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
''' morphology_analysis.py  -   compare the habits of two crystal structures.

    Determines the aspect ratios of the BFDH morphologies of the two crystals
    and plots them on a chart, showing which class of habit they are likely to
    possess.
'''

import argparse
import os
import matplotlib
matplotlib.use('Agg')
from matplotlib import pyplot as plt

from ccdc.descriptors import CrystalDescriptors
from ccdc.io import EntryReader

arg_parser = argparse.ArgumentParser(description=__doc__)
arg_parser.add_argument(
    'refcodes', nargs='+', help='Refcodes of structures to analyse'
)
arg_parser.add_argument(
    '-o', '--output', default=os.path.join('.', 'morphology_aspects.png'), help='Output chart PNG filename.'
)
args = arg_parser.parse_args()

#  Set up plot
fig = plt.figure(figsize=(10., 10.))
ax = plt.axes()
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.set_xlabel('M:L Ratio', fontsize=16)
ax.set_ylabel('S:M Ratio', fontsize=16)
ax.annotate('Needle', (0.23, 0.95), fontsize=16)
ax.annotate('Block', (0.73, 0.95), fontsize=16)
ax.annotate('Lath', (0.23, 0.45), fontsize=16)
ax.annotate('Plate', (0.73, 0.45), fontsize=16)
plt.plot([0.5, 0.5], [0, 1], c='black', linewidth=2)
plt.plot([0, 1], [0.5, 0.5], c='black', linewidth=2)
ptc = fig.patch
ptc.set_facecolor('white')


def make_plot(asp, colour, refcode):
    #  This function plots the ratios
    y, x = get_ratios(asp)
    plt.scatter(x, y, c=colour, s=40)
    ax.annotate('{}'.format(refcode), (x + 0.007, y - 0.007), fontsize=10)


def get_ratios(aspect):
    #  This function returns the side length ratios small:medium and medium:large
    s_m = aspect[0] / aspect[1]  # Small to medium side lengths ratio
    m_l = aspect[1] / aspect[2]  # Medium to long side lengths ratio
    return s_m, m_l


def get_form(aspect):
    #  In this function we designate a shape from the aspect ratio
    s_m, m_l = get_ratios(aspect)
    if s_m < 0.5:
        if m_l < 0.5:
            return 'Lath'
        else:
            return 'Plate'
    elif m_l < 0.5:
        return 'Needle'
    else:
        return 'Block'


if __name__ == '__main__':
    # Read structures from the CSD
    reader = EntryReader('CSD')
    if not os.path.exists(os.path.dirname(args.output)):
        os.makedirs(os.path.dirname(args.output))

    # Read in two crystal polymorphs of the drug Nabumetone
    for refcode in args.refcodes:
        crystal = reader.crystal(refcode)
        morphology = CrystalDescriptors.Morphology(crystal)
        bounding_box = morphology.oriented_bounding_box
        #  By calculating the bounging box of the morphologies, we can get the aspect ratios
        aspect = [
            bounding_box.minor_length, bounding_box.median_length, bounding_box.major_length
        ]
        # The aspect ratios can be used to designate a shape for the crystal
        shape = get_form(aspect)
        print('{}: {}'.format(refcode, shape))
        # This can be plotted
        make_plot(aspect, 'red', refcode)
        # Using the plotting from Morphology.py we can also overlay the morphologies in 3D
        plt.savefig(args.output)

Charting crystal morphologies

This example draws an image of the morphology and tables of plane angles and of face relative areas of morphologies of one or more crystals of the CSD. An example for HXACAN is in ../images/hxacan_morphology.html.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
'''morphology_plot.py   -   make an image of a morphology.
'''

import argparse
import codecs
import itertools
import io as sys_io
import os
import matplotlib
matplotlib.use('Agg')
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

from ccdc.io import EntryReader
from ccdc.descriptors import CrystalDescriptors
from ccdc.descriptors import GeometricDescriptors

# Visualisation functions


def wireframe(face, ax):
    #  Plot wireframe of morphology in matplotlib
    for n, edge in enumerate(face.edges):
        Axes3D.plot(ax,
                    [coord[0] for coord in edge],
                    [coord[1] for coord in edge],
                    [coord[2] for coord in edge],
                    c='black',
                    linewidth=1.5)


def skin(face, colour, ax):
    #  Mesh the facets for matplotlib
    for edge in face.edges:
        vertices = [(edge[0], edge[1], face.centre_of_geometry)]
        Axes3D.add_collection3d(ax, Poly3DCollection(
            vertices, color=colour, linewidth=0, alpha=0.1))


def add_labels(morphology, face, ax):
    #  Add labels to facets
    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')


def generate_plot(morphology, colour):
    #  Plot the morphology in matplotlib
    # Set up Axes
    fig = plt.figure(figsize=(10., 10.))  # 3D graph instance
    ax = fig.add_subplot(111, projection='3d')  # 3D Axes
    ax.set_xlim3d(-20, 20)
    ax.set_ylim3d(-20, 20)
    ax.set_zlim3d(-20, 20)
    ax.set_xlabel('X axis')
    ax.set_ylabel('Y axis')
    ax.set_zlabel('Z axis')
    ax.grid(False)
    plt.axis('off')

    for facet in morphology.facets:
        wireframe(facet, ax)
        skin(facet, colour, ax)
        add_labels(morphology, facet, ax)
    return ax, plt


# Analysis functions


def compare_angles(morphology, format='html'):
    # Print a table of every combination of facets and their dihedral angle
    indices = {}
    for f in morphology.facets:
        indices[f.miller_indices.hkl] = GeometricDescriptors.Plane.from_points(*f.coordinates)
    # Every pair of facets to compare
    combinations = list(itertools.combinations([k for k in indices], 2))
    if format == 'html':
        return '<table><tr><th>Plane1</th><th>Plane2</th></th><th>Angle</th>%s</table>' % (
            '\n'.join('<tr><td>%s</td><td>%s</td><td align="right">%.2f</td></tr>' %
                      (c[0], c[1], 180 - (indices[c[0]].plane_angle(indices[c[1]])))
                      for c in combinations)
        )
    for combination in combinations:
        angle = 180 - (indices[combination[0]]).plane_angle(indices[combination[1]])
        print('planes {} & {}: angle = {}'.format(combination[0], combination[1], angle))


def print_areas(morphology, format='html'):
    '''Print a table of every face and its relative area.'''
    if format == 'html':
        return '<table><tr><th>Indices</th><th>Area</th></tr>%s</table>' % (
            '\n'.join('<tr><td>%s</td><td align="right">%.2f</td></tr>' %
                      (face.miller_indices.hkl, morphology.relative_area(face.miller_indices))
                      for face in morphology.facets
                      )
        )
    for face in morphology.facets:
        print('{}: relative area = {}'.format(
            face.miller_indices.hkl, round(morphology.relative_area(face.miller_indices), 3)))


def run_all(morphology):
    # Run all of the functions in this script
    ax, plt = generate_plot(morphology, (0, 0, 1, 0.1))
    compare_angles(morphology)
    print_areas(morphology)
    return ax, plt


if __name__ == '__main__':
    arg_parser = argparse.ArgumentParser(description=__doc__)
    arg_parser.add_argument(
        'refcodes', nargs='+', help='refcodes of CSD structures to plot'
    )
    arg_parser.add_argument(
        '-o', '--output', default=os.path.join('.', 'morphologies.html'),
        help='HTML file for the output report.'
    )

    args = arg_parser.parse_args()
    reader = EntryReader('CSD')  # Read structures from the CSD
    if not os.path.exists(os.path.dirname(args.output)):
        os.makedirs(os.path.dirname(args.output))

    lines = [
        '<style>',
        '.box { margin:0; }',
        '.left { float: left; width: 50%; overflow: hidden; }',
        '.middle { float: left; width: 25%; overflow: hidden; }',
        '.right { float: left; width: 20%; overflow: hidden; }',
        '</style>',
        '']

    for refcode in args.refcodes:
        cryst = reader.crystal(refcode)  # Open the crystal from CSD entry SUCACB03
        morphology = CrystalDescriptors.Morphology(cryst)  # Create a morphology object
        ax, plt = generate_plot(morphology, (0, 0, 1, 0.1))
        svg = sys_io.BytesIO()
        plt.savefig(svg, format='svg')
        lines.append('<div class="box">')
        lines.append('<div class="left">%s</div>' % svg.getvalue().decode('utf-8'))
        lines.append('<div class="middle">%s</div>' % compare_angles(morphology, format='html'))
        lines.append('<div class="right">%s</div>' % print_areas(morphology, format='html'))
        lines.append('</div>')

        with codecs.open(args.output, 'w', encoding='utf-8') as writer:
            for line in lines:
                writer.write('%s%s' % (line, os.linesep))

Functional groups of a crystal’s HBonds

This writes a csv file containing the hydrogen bonding atoms of crystals with the functional groups of the donor and acceptor atoms.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
'''
    functional_groups_of_hbonds.py  -   produce a spreadsheet of functional groups of hbonds
'''
################################################################################

import sys
import csv
import argparse

from ccdc.io import CrystalReader
from ccdc.descriptors import CrystalDescriptors

################################################################################

class Runner(argparse.ArgumentParser):
    '''Defines and parses arguments, runs the job.'''
    def __init__(self):
        super(self.__class__, self).__init__(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
        self.add_argument(
            'gcd_file',
            help='GCD file of crystal structures to analyse'
        )
        self.add_argument(
            '-o', '--output-file', '--output_file', default=None,
            help='Output csv file [stdout]'
        )
        self.args = self.parse_args()

    def run(self):
        '''Runs the job.'''
        calculator = CrystalDescriptors.HBondCoordination()
        reader = CrystalReader(self.args.gcd_file)
        if self.args.output_file:
            f = open(self.args.output_file, 'w')
        else:
            f = sys.stdout
        headers = [
            'Identifier', 'Donor residue', 'Donor functional group', 'Donor atom',
            'Acceptor residue', 'Acceptor functional group', 'Acceptor atom'
        ]
        writer = csv.writer(f)
        writer.writerow(headers)
        def get_residue(crystal, label):
            for i, c in enumerate(crystal.molecule.components):
                try:
                    c.atom(label)
                    return 'RES%d' % (i+1)
                except RuntimeError:
                    pass
            raise RuntimeError('Ambiguous atom label %s' % label)


        for c in reader:
            if c.has_disorder:
                continue
            predictions = calculator.predict(c)
            for hb in c.hbonds():
                if hb.atoms[1] in hb.atoms[0].neighbours:
                    don_at, acc_at = hb.atoms[0], hb.atoms[-1]
                else:
                    don_at, acc_at = hb.atoms[-1], hb.atoms[0]
                don, acc = predictions.functional_groups_of_hbond(hb)
                if don is None:
                    don = '%s Nothing Unidentified group Nothing' % don_at.label
                if acc is None:
                    acc = '%s Nothing Unidentified group Nothing' % acc_at.label
                don_bits = don.split()
                acc_bits = acc.split()
                row = [
                    c.identifier,
                    get_residue(c, don_bits[0]),
                    ' '.join(don_bits[2:-1]),
                    don_bits[0],
                    get_residue(c, acc_bits[0]),
                    ' '.join(acc_bits[2:-1]),
                    acc_bits[0]
                ]
                writer.writerow(row)

################################################################################

if __name__ == '__main__':
    Runner().run()

A simple spreadsheet of hbond coordination likelihoods

This example iterates over a GCD file (or any other crystal file) creating a csv file of hbond coordination likelihoods.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
'''
    hbond_coordination_prediction_spreadsheet.py - calculate hbond coordination predictions for a set of crystals

    This makes a spreadsheet of hbond coordination predictions containing columns indicating predicted probability and whether or not that number of hbonds was observed.
'''
#############################################################################################

import argparse

from ccdc.io import CrystalReader
from ccdc.descriptors import CrystalDescriptors
from ccdc.utilities import CSVWriter

#############################################################################################

class Runner(argparse.ArgumentParser):
    '''Defines arguments, parses arguments and runs the analysis.'''
    def __init__(self):
        super(self.__class__, self).__init__(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
        self.add_argument(
            'gcd_file', type=str,
            help='GCD file of crystal structures to analyse'
        )
        self.add_argument(
            '-o', '--output-file', '--output_file', default=None,
            help='Output csv file [stdout]'
        )
        self.add_argument(
            '-m', '--model-path', '--model_path', default=None,
            help='Path to a directory of coordination models'
        )
        self.args = self.parse_args()

    def run(self):
        settings = CrystalDescriptors.HBondCoordination.Settings()
        if self.args.model_path is not None:
            settings.coordination_models_path = self.args.model_path
        calculator = CrystalDescriptors.HBondCoordination(settings=settings)
        reader = CrystalReader(self.args.gcd_file)

        outfile = self.args.output_file if self.args.output_file else None

        headers = ['Identifier', 'Molecule', 'Atom', 'Functional group', 'D/A']
        for i in range(7):
            headers.extend(['%s %d likelihood' % ('>=' if i == 6 else '=', i), 'Observed'])

        with CSVWriter(file_name=outfile, header=headers) as writer:
            for c in reader:
                if c.has_disorder:
                    continue
                predictions = calculator.predict(c)
                for r in predictions.to_csv(separator=None):
                    writer.write_row(r)

#############################################################################################

if __name__ == '__main__':
    Runner().run()

Exploring hydrogen bond coordination models

This example builds a coordination model as a tree of possible coordinations for the donors and acceptors of a crystal. This is only one of many possible models that can be constructed from the API.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
"""
    hbond_coordination_tree.py  -   create a tree of possible hbond coordinations with attached probabilities.

    The tree enumerates each possible coordination for each donor atom. Subtrees from this enumerates each possible coordination for each acceptor atom consistent with the donations previously defined.  For example in the subtree where the first donor donates once and the second donor donates once, the sum of the acceptor atoms' coordination must sum to two.

    For convenience, the created tree will be returned, so invoking this interactively allows one to explore the tree data structure.
"""

import argparse
import os

from ccdc.descriptors import CrystalDescriptors
from ccdc import io
from ccdc.diagram import DiagramGenerator
from ccdc.utilities import CSVWriter, HTMLReport, html_table

csd = io.EntryReader('csd')


class Node(object):
    """A node in a hydrogen bond propensity tree.

    :param label: (:obj:`str`) A label for the node
    :param is_donor: (:obj:`boolean`) Whether the node is a donor.
    :param coordination: (:obj:`float`) This node's coordination score.
    :param probability: (:obj:`float`) The hydrogen bond probability of this node.
    :param cumulative_probability: (:obj:`float`)  The cumulative hydrogen bond probability of this
        node and all the nodes on the path to it.
    :param is_observed: (:obj:`boolean`) Whether the node is on the observed coordination path.
    :param current_donations: (:obj:`list`) The nodes which are donors for this node.
    :param remaining_donors: (:obj:`list`) The nodes which are donors for this node but have not
        yet been factored into the cumulative probability.
    :param remaining_acceptors: (:obj:`list`) Remaining acceptors for this node.
    :param predictions: (:obj:`list`) The probability predictions for this node.
    """

    def __init__(self, label, is_donor, coordination, probability, cumulative_probability,
                 is_observed, current_donations, remaining_donors, remaining_acceptors, predictions):
        """Holds all the info for a decision tree node."""
        self.label = label
        self.is_donor = is_donor
        self.coordination = coordination
        self.probability = probability
        self.cumulative_probability = cumulative_probability * self.probability
        self.is_observed = is_observed
        self.remaining_donors = remaining_donors
        self.remaining_acceptors = remaining_acceptors
        self.current_donations = current_donations
        self.predictions = predictions
        self.children = []
        self._expand()

    def _expand(self):
        """Expands child nodes."""
        if self.remaining_donors:
            preds = self.predictions.predictions_for_label(self.remaining_donors[0])
            self.children = [
                Node(self.remaining_donors[0], True, p, preds[1][p], self.cumulative_probability,
                     p == preds[0], p + self.current_donations,
                     self.remaining_donors[1:], self.remaining_acceptors, self.predictions)
                for p in preds[1]
                if preds[1][p] != 0.0 and preds[1][p]*self.cumulative_probability > 0.0
            ]
        elif self.remaining_acceptors:
            preds = self.predictions.predictions_for_label(self.remaining_acceptors[0], 'acceptor')
            self.children = [
                Node(self.remaining_acceptors[0], False, p, preds[1][p],
                     self.cumulative_probability, p == preds[0], self.current_donations - p,
                     self.remaining_donors, self.remaining_acceptors[1:], self.predictions)
                for p in preds[1]
                if preds[1][p] != 0.0 and p <= self.current_donations
                and preds[1][p] != 0.0 and preds[1][p] * self.cumulative_probability > 0.0
            ]

    @property
    def is_leaf(self):
        """:obj:`boolean`: Whether or not this node is a leaf node."""
        return self.children == []

    @property
    def observed_path(self):
        """:obj:`list` of :obj:`Node`: All nodes corresponding to the crystal's observed coordination."""
        if self.is_leaf:
            if self.is_observed:
                return [self]
        else:
            for k in self.children:
                if k.is_observed:
                    return [self] + k.observed_path
        print(self.label, self.is_observed, self.coordination)
        print([(k.is_observed, k.label, k.coordination) for k in self.children])
        raise RuntimeError('No observed child')

    @property
    def leaves(self):
        """:obj:`list` of :obj:`Node`: All the leaf nodes in the tree."""
        if self.is_leaf:
            yield self
        else:
            for k in self.children:
                for l in k.leaves:
                    yield l

    def walk(self):
        """Walks over the tree yielding all nodes.

        :returns: (:obj:`Node`) Yields all the nodes in the tree.
        """
        yield self
        for k in self.children:
            for x in k.walk():
                yield x

    def _set_depth(self, d):
        """Recursively sets the depth of the tree.

        :param d: (:obj:`int`) The current depth of the tree.
        """
        self.depth = d
        for k in self.children:
            k._set_depth(d+1)

    def _prune(self):
        """Remove paths where the leaf has still got unaccounted for donations."""
        for k in self.children:
            k._prune()
        to_go = [(i, k) for i, k in enumerate(self.children) if k.is_leaf and k.current_donations]
        to_go.sort(reverse=True)
        for i, k in to_go:
            self.children.pop(i)

    def _accumulate_das(self, dons, accs):
        """Sets in the node the accumulated coordinations.

        This will be attributes, dons and accs, contaiming labels of the donor or acceptor with a
        multiplicity given by the coordination.

        :param dons: (:obj:`list`) A list of donor atoms.
        :param accs: (:obj:`list`) A list of acceptor atoms.
        """
        if self.is_donor:
            dons = dons + [self.label] * self.coordination
        else:
            accs = accs + [self.label] * self.coordination
        if self.is_leaf:
            self.donors = dons
            self.acceptors = accs
        else:
            for k in self.children:
                k._accumulate_das(dons, accs)

    def _set_cumulative_prob(self, prob):
        """There are various different ways of setting the cumulative probability.
        Here we have chosen to ignore the probability associated with coordination scores of 0.

        :param prob: (:obj:`float`) The probability to factor in to the total.
        """
        # Ignore 0-coordinated acceptors
        if not self.is_donor and self.coordination == 0:
                self.cumulative_probability = prob
        else:
            self.cumulative_probability = self.probability * prob
        for k in self.children:
            k._set_cumulative_prob(self.cumulative_probability)

    def _set_full_cumulative_prob(self, prob):
        """Sets the full cumulative probability including the 0-coordination probabilities.

        :param prob: (:obj:`float`) The probability to factor in to the total.
        """
        self.full_cumulative_probability = self.probability * prob
        for k in self.children:
            k._set_full_cumulative_prob(self.full_cumulative_probability)

    def _accumulate_coordinations(self, accumulated):
        """Accumulates the coordination in a list.

        This is to ensure leaf nodes have the complete list of coordinations contributing.

        :param accumulated: (:obj:`list`) The coordinations to add to the total.
        """
        self.accumulated_coordination = accumulated + [self.coordination]
        for k in self.children:
            k._accumulate_coordinations(self.accumulated_coordination)


class RootNode(Node):
    """The root node of the tree.

    :param donors: (:obj:`list`) The donor atoms.
    :param acceptors: (:obj:`list`) The acceptor atoms.
    :param predictions: (:obj:`list`) The probability predictions.
    :param children: (:obj:`Node`) The child nodes of the root.
    """

    def __init__(self, donors, acceptors, predictions, children):
        self.label = 'ROOT'
        self.donors = donors
        self.acceptors = acceptors
        self.is_donor = None
        self.coordination = 0
        self.probability = 1.0
        self.cumulative_probability = 1.0
        self.is_observed = True
        self.remaining_donors = donors
        self.remaining_acceptors = acceptors
        self.current_donations = 0
        self.predictions = predictions
        self.children = children

        for k in self.children:
            k._expand()

        self._prune()
        self._set_depth(0)
        self._accumulate_das([], [])
        self._set_cumulative_prob(1.0)
        self._set_full_cumulative_prob(1.0)

        for k in self.children:
            k._accumulate_coordinations([])

    @staticmethod
    def make_tree(crystal):
        """Makes the tree from the crystal.

        :param crystal: (:obj:`ccdc.crystal.Crystal`) A CSD Python API Crystal object.
        :returns: (:obj:`RootNode`) A Hydrogen Bond Propensity tree for the crystal.
        """
        coordination_calc = CrystalDescriptors.HBondCoordination()
        predicts = coordination_calc.predict(crystal)
        dons = [d.label for d in coordination_calc.donor_atoms]
        accs = [a.label for a in coordination_calc.acceptor_atoms]
        preds = predicts.predictions_for_label(dons[0])
        return RootNode(dons, accs, predicts, [
            Node(dons[0], True, p, preds[1][p], 1.0, p == preds[0], p, dons[1:], accs, predicts)
            for p in preds[1] if preds[1][p] != 0.0
        ])


class Runner(argparse.ArgumentParser):
    """Handle command-line arguments and generate the reports."""

    def __init__(self):
        """Set up the argument parser."""
        super(self.__class__, self).__init__(description=__doc__,
                                             formatter_class=argparse.ArgumentDefaultsHelpFormatter)
        self.add_argument('refcodes', nargs='+', help='CSD Refcodes or files')
        self.add_argument('-o', '--output',
                          help='The output file to use. HTML or CSV files can be written.')
        self.args = self.parse_args()


        self.output_format = None
        """:obj:`str` An identifier for the output format to write."""

        self.output_file = None
        """:obj:`str` The full file name of the output file."""

        self.output_dir = None
        """:obj:`str` The directory to hold output files."""

        self.diagram_generator = None
        """:obj:`ccdc.diagram.DiagramGenerator` A diagram generator."""

        if not self.args.output:
            self.output_format = 'stdout'
        else:
            if not (self.args.output.endswith('html') or self.args.output.endswith('csv')):
                raise Exception('Cannot recognise format for output file %s! .html or .csv are valid.'
                                % self.args.output)

            self.output_file = os.path.abspath(self.args.output)
            self.output_dir = os.path.dirname(self.output_file)

            # Make sure the output directory exists
            if not os.path.exists(self.output_dir):
                os.makedirs(self.output_dir)

            if self.output_file:
                self.diagram_generator = DiagramGenerator()
                self.diagram_generator.settings.return_type = 'SVG'

            self.output = None
            if self.args.output.endswith('html'):
                self.output_format = 'html'
            elif self.args.output.endswith('csv'):
                self.output_format = 'csv'

    def run(self):
        """Run the reporting on all specified crystals.

        :returns: (:obj:`RootNode`) A Hydrogen Bond Propensity tree for the last crystal processed.
        """
        for refcode_or_file in self.args.refcodes:
            if os.path.isfile(refcode_or_file):
                for crystal in io.CrystalReader(refcode_or_file):
                    tree = self.run_one(crystal)
            else:
                tree = self.run_one(csd.crystal(refcode_or_file))

        # Write the HTML footer for HTML reports
        if self.output_format == 'html':
            with self.get_output_writer() as report:
                report.write_footer()

        return tree

    def get_output_writer(self, identifier=''):
        """Get an output writer for the given identifer.

        This will reopen and close the files after each iteration so problems any structure
        don't result in all the data so far being lost.

        :param identifier: (:obj:`str`) The identifier of the structure to write for.
        :returns: (:obj:`file`) A context handler file-like object that can write output,
            appropriate to the output format specified.
        """
        if self.output_format == 'html':
            return HTMLReport(file_name=self.output_file,
                              report_title='Hydrogen Bond Propensity Report',
                              append=True)
        elif self.output_format == 'csv':
            out_path = os.path.dirname(self.output_file)
            basename, ext = os.path.splitext(os.path.basename(self.output_file))
            output_file_path = os.path.join(out_path, '%s_%s%s' % (basename, identifier, ext))
            return CSVWriter(file_name=output_file_path)
        elif self.output_format == 'stdout':
            return CSVWriter(stdout=True)

    def run_one(self, crystal):
        """Calculates the tree for one crystal and writes results to output.

        :param crystal: (:obj:`ccdc.crystal.Crystal`) A CSD Python API Crystal object.
        :returns: (:obj:`RootNode`) A Hydrogen Bond Propensity tree.
        """
        tree = RootNode.make_tree(crystal)
        if self.output_format == 'html':
            self.write_html(tree, crystal)
        else:
            self.write_csv(tree, crystal)

        return tree

    @staticmethod
    def csv_headers(tree):
        """Generate CSV headers for the given tree.

        :param tree: (:obj:`RootNode`) A Hydrogen Bond Propensity tree.
        :returns: (:obj:`list`) A list of CSV column headers for the given tree.
        """
        headers = ['Identifier']
        headers += ['D %s' % tree.donors[i] for i in range(len(tree.donors))]
        headers += ['A %s' % tree.acceptors[i] for i in range(len(tree.acceptors))]
        headers += ['Probability', 'Full Probability', 'Observed']
        return headers

    def write_csv(self, tree, crystal):
        """Write tree data to CSV.

        :param tree: (:obj:`RootNode`) A Hydrogen Bond Propensity tree.
        :param crystal: (:obj:`ccdc.crystal.Crystal`) A CSD Python API Crystal object.
        """
        observed_path = tree.observed_path
        leaves = [(n.cumulative_probability, n) for n in tree.leaves]
        leaves.sort(reverse=True)

        with self.get_output_writer(crystal.identifier) as writer:
            writer.write_header(self.csv_headers(tree))
            for probability, leaf in leaves:
                row = [crystal.identifier]
                row += list(leaf.accumulated_coordination)
                row += [probability, leaf.full_cumulative_probability, leaf == observed_path[-1]]
                writer.write_row(row)

    @staticmethod
    def propensity_data(leaves, observed_path):
        """Generate rows of HTML output.

        :param leaves: (:obj:`list`) A list of (:obj:`float`, :obj:`Node`) tuples of hydrogen bond
            probabilities and the leaves of the probability tree they correspond to.
        :param observed_path: (:obj:`list`) A list of :obj:`Node`s representing the observed
            coordination of the crystal.
        :returns: (:obj:`list`) A list of the propensity data for each leaf passed in.
        """
        for leaf_data in leaves:
            leaf = leaf_data[1]
            probability = leaf_data[0]
            yield [str(i) for i in leaf.accumulated_coordination] +\
                  ['%.4f' % probability, '%.4f' % leaf.full_cumulative_probability,
                   (leaf == observed_path[-1])]

    def write_html(self, tree, crystal):
        """Writes the given crystal's report to HTML.

        :param tree: (:obj:`RootNode`) A Hydrogen Bond Propensity tree.
        :param crystal: (:obj:`ccdc.crystal.Crystal`) A CSD Python API Crystal object.
        """
        # Sort the data by reverse probability
        leaves = [(n.cumulative_probability, n) for n in tree.leaves]
        leaves.sort(reverse=True)

        # Generate a diagram
        diagram = self.diagram_generator.image(crystal).replace(
            'Qt SVG Document', 'Diagram for %s' % crystal.identifier)

        with self.get_output_writer() as report:
            report.write_section_header(headline=crystal.identifier)
            report.write_figure_data(diagram, caption='Diagram for %s' % crystal.identifier,
                                     figure_id='diagram_%s' % crystal.identifier)
            table_header = ['D %s' % x for x in tree.donors] +\
                           ['A %s' % a for a in tree.acceptors] +\
                           ['Probability', 'Full Probability', 'Observed']
            report.write(html_table(data=self.propensity_data(leaves,
                                                              observed_path=tree.observed_path),
                                    header=table_header,
                                    table_id='data_%s' % crystal.identifier,
                                    caption='Propensity data for %s' % crystal.identifier))


if __name__ == '__main__':
    t = Runner().run()

Charting hydrogen bond propensities

This example creates a chart of hydrogen bond propensity groupings, plotting hydrogen bond score against hbond coordination score. Where there are structures to the right of and lower than the observed structure, there is a greater chance of the existence of polymorphic structures, for example:

YIGPIO03 hydrogen bond propensity chart

YIGPIO03 Hydrogen Bond Propensities

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at http://www.ccdc.cam.ac.uk/Community/Pages/Licences/v2.aspx
#
# This permission notice and the following statement of attribution must be
# included in all copies or substantial portions of this script.
#
'''
    hydrogen_bond_propensities.py - performs an hbp calculation and produces a chart of propensity vs coordination score
'''
################################################################################

import argparse
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

import ccdc.io
from ccdc.descriptors import CrystalDescriptors

################################################################################


def run(refcode):
    '''Runs the job.'''
    hbp = CrystalDescriptors.HBondPropensities()
    reader = ccdc.io.CrystalReader('csd')
    crystal = reader.crystal(refcode)
    print('Hydrogen Bond Propensity calculation for %s\n' % refcode)
    hbp.set_target(crystal)
    hbp.match_fitting_data(count=300)
    hbp.analyse_fitting_data()
    hbp.perform_regression()
    propensities = hbp.calculate_propensities()
    groups = hbp.generate_hbond_groupings()
    observed_groups = hbp.target_hbond_grouping()
    make_chart(groups, observed_groups, crystal.identifier)


def make_chart(groups, observed_groups, identifier):
    '''Makes a chart of the hbond groupings, plotting hbond score against coordination score.'''
    figure = plt.scatter([group.hbond_score for group in groups],
                         [group.coordination_score for group in groups],
                         s=50, label='Putative Structures')
    ax2 = plt.scatter(observed_groups.hbond_score,
                      observed_groups.coordination_score,
                      c='red', marker='*', s=250, label='Observed Structure')
    plt.origin = 'upper'
    plt.xlim(0, 1.0)
    plt.ylim(-1.0, 0)
    plt.legend(scatterpoints=1)
    ax = plt.gca()
    ax.xaxis.tick_top()
    ax.yaxis.tick_left()
    ax.set_xlabel('Mean H-Bond Propensity')
    ax.set_ylabel('Mean H-Bond Coordination')
    ax.xaxis.set_label_position('top')
    fig_name = '%s_standard_chart' % identifier
    plt.savefig('%s.png' % fig_name, format='png')
    plt.savefig('%s.svg' % fig_name, format='svg')


################################################################################

if __name__ == '__main__':
    parser = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter, description=__doc__)
    parser.add_argument('refcode', help='Refcode of structure to be analysed', type=str)

    args = parser.parse_args()

    run(args.refcode)