Particle examples

Note

The particle API is available only to CSD-Particle users.

Slip planes

This example shows how to extract information about possible slip planes and save the data to a csv file.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at https://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.
#
#
# 2020-09-04: made available by the Cambridge Crystallographic Data Centre
#

'''
slip_planes.py - extract information about possible slip planes and can save it to csv.

The script uses the CCDC particle module to extract information, and,
optionally save the data to a csv file and select the desired separation threshold

Example:
    python3 slip_planes.py HXACAN

'''

###########################################################################
import argparse
import csv
from itertools import zip_longest
import os

from ccdc.io import CrystalReader
from ccdc import particle


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

def parse_arguments():
    parser = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.RawTextHelpFormatter)
    parser.add_argument(
        'refcode',
        help='The REFCDOE of a CSD entry, eg. HXACAN28.'
    )
    parser.add_argument(
        '-s', '--separation',
        type=float,
        default=0.0,
        help='Set desired separation threshold, default is set to 0.'
    )
    parser.add_argument(
        '-csv', '--csv_output',
        action='store_true',
        help='Save output as csv.'
    )
    return parser.parse_args()


def calculate_slip_planes(crystal, threshold):
    """
    Calculates slip planes and set separation threshold

    :param crystal: The crystal whose slip planes we want to calculate
    :type crystal: class:`ccdc.crystal.Crystal`
    :param threshold: values between [-2 , 0] to select separation between planes
    :type threshold: float
    :return: planes object
    """
    print("Calculating possible slip planes...")
    planes = particle.SlipPlanes(crystal)
    planes.set_separation_threshold(threshold)

    return planes


def print_plane_info(plane, count):
    """
    Prints out information about the given plane

    :param plane: Plane of interest
    :type plane: class: `ccdc_addop.SlipPlanes.SlipPlane` object
    :param count: number in planes list
    :type count: interger
    """
    print(f"\n  Plane number {count}")
    print(f"    Slip plane orientation     : {plane.orientation.hkl}")
    print(f"    Offset                     : {round(plane.offset,3)}A")
    print(f"    Slab Separation            : {round(plane.slab_separation,3)}A")
    print(f"    Repeat Distance            : {round(plane.repeat_distance,3)}A")
    print(f"    Perpendicular planes       : {[p_plane.orientation.hkl for p_plane in plane.perpendicular_slip_planes]}")
    print(f"    Hydrogen bond across plane : {plane.h_bonded}")


def save_to_csv(planes_dict, refcode):
    """
    Outputs the dictionary as a csv to the current directory
    :param planes_dict: Dictionary of lists containing all the slip planes
    :type plane_dict: dict
    :param refcode: refcode of the structure being saved
    :type refcode: str
    """
    transposed_data = list(zip_longest(*planes_dict.values()))
    filename = f"{refcode}_SlipPlanes.csv"
    try:
        with open(filename, 'w', newline='') as f:
            writer = csv.writer(f)
            writer.writerow(planes_dict.keys())
            for items in transposed_data:
                writer.writerow(items)
        print(f"`\n File saved in : {os.getcwd()}\\{filename}")
    except IOError as e:
        print(str(e))
        print("Can't write to this folder")


def main(args):
    """ Driving function"""
    crystal = CrystalReader('CSD').crystal(args.refcode)
    print(f"Crystal Selected - {args.refcode}")
    planes = calculate_slip_planes(crystal, args.separation)

    divider = "######" * 5
    print(f"\n{divider} Slip Planes {divider}")
    print(f"Number of slip planes with a separation > {planes.separation_threshold()} : {len(planes)}")
    if planes.hbond_dimensionality < 0:
        print(f"No hydrogen bonds")
    else:
        print(f"Hydrogen bond dimensionality of crystal : {planes.hbond_dimensionality}-D")

    if len(planes) > 0:
        print(f"Possible slip planes : {[plane.orientation.hkl for plane in planes]}")
        print(f"\n {divider} Individual plane analysis Details {divider}")
        for i, plane in enumerate(planes):
            print_plane_info(plane, i)

    else:
        print(f"No planes found with a separation > {args.separation}A")

    print(f"\n {divider} Done {divider}")

    if args.csv_output:
        save_to_csv(planes.as_dict(), args.refcode)




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

if __name__ == '__main__':
    args = parse_arguments()
    main(args)

Surface

This example shows how to extract information and visualize the surface of a crystal plane.

#!/usr/bin/env python
#
# This script can be used for any purpose without limitation subject to the
# conditions at https://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.
#
# 2020-06-11: created by Alex Moldovan, the Cambridge Crystallographic Data Centre
# 2020-06-16: made available by the Cambridge Crystallographic Data Centre
#

'''
surface.py - extract information and plot contour of the surface.

The script uses the CCDC particle module to extract information, and,
optionally, visualize the topology, of a crystal plane surface.

Example:
    python surface.py NALCYS23 (101)

'''
###########################################################################

import argparse
import os
import re
from typing import Tuple, Dict

try:
    import numpy as np
    from scipy.spatial.transform import Rotation
    from scipy.interpolate import griddata
    from scipy.ndimage import gaussian_filter
    import matplotlib.pyplot as plt
except ImportError as exc:
    print(str(exc))
    print('This is needed to visualize the surface topology\n')

from ccdc.io import CrystalReader
from ccdc import particle


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


def m_indices(input_string: str) -> Tuple:
    """To convert the hkl argument into tuple used to create Surface object.

    Example:
        input: "(110)"
        output: tuple([1,1,0])
    """
    mi_regex = r'\((-?\d)(-?\d)(-?\d)\)'
    mm = re.match(mi_regex, input_string)
    if mm is None:
        raise argparse.ArgumentTypeError('"{}" is not of the form "(hkl)"' \
                                         .format(input_string))
    return tuple(map(int, mm.group(1, 2, 3)))


def parse_arguments():
    parser = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.RawTextHelpFormatter)
    parser.add_argument(
        'refcode',
        help='The Refcode of a CSD entry, eg. NALCYS23.'
    )
    parser.add_argument(
        'hkl',
        type=m_indices,
        help='The Miller indices of the plane, eg. (101).'
    )
    parser.add_argument(
        '-p', '--plot_topology',
        action='store_true',
        help='Plot topology as contour map.'
    )
    parser.add_argument(
        '-n', '--node_data',
        action='store_true',
        help='Save node data as output file in directory.'
    )
    return parser.parse_args()


def _rotate_to_xy(normal: np.ndarray, cleave_edges: list) -> np.ndarray:
    """
    Rotates the plane to fit xy plane

    :param normal: numpy.array vector normal of the plane
    :param cleave_edges: list of interception edges across BFDH
    :return: numpy.array rotated edges to given alignment axis Default = [0,0,1]
    """
    alignment_axis = np.array([0, 0, 1])
    if np.array_equal(normal, alignment_axis):
        return np.array(cleave_edges)

    cross = np.cross(normal, alignment_axis)
    norm = np.linalg.norm(cross)
    if norm:
        cross = cross / norm
    theta = np.arccos((np.dot(normal, alignment_axis)))
    rot = Rotation.from_rotvec(theta * cross)
    return np.array([rot.apply(x) for x in cleave_edges])


def _generate_interpolated_grid_data(nodes: np.ndarray, cell_limits: Dict, method: str = 'linear',
                                     grid_spacing: float = 0.2) -> np.ndarray:
    """
    Uses nodes to generate an interpolated grid

    :param list nodes: coordinates of node positions
    :param str method: default 'linear'. Look up scipy.interpolate.griddata for more options.
    :return:
    """
    if not isinstance(nodes, np.ndarray):
        nodes = np.array(nodes)

    x_line = np.arange(cell_limits['x_min'], cell_limits['x_max'], grid_spacing)
    y_line = np.arange(cell_limits['y_min'], cell_limits['y_max'], grid_spacing)

    x, y = np.meshgrid(x_line, y_line)

    coord_xy = [[nodes[:, 0][i], nodes[:, 1][i]] for i in range(len(nodes[:, 0]))]
    if len(nodes) < 10:
        method = 'nearest'
    grid = griddata(np.array(coord_xy), nodes[:, 2], (x, y), method=method)

    return grid


def _node_limits(nodes: np.ndarray) -> Dict:
    """ Calculates the min and max of all nodes """

    if not isinstance(nodes, np.ndarray) and isinstance(nodes, list):
        nodes = np.array(nodes)

    x_min = nodes[:, 0].min()
    y_min = nodes[:, 1].min()
    x_max = nodes[:, 0].max()
    y_max = nodes[:, 1].max()

    return {'x_min': x_min, 'y_min': y_min, 'x_max': x_max, 'y_max': y_max}


def plot_contour(array: np.ndarray, refcode: str, hkl: str) -> None:
    """ Plots a contour map based on the nodes of the surface topology """
    fig = plt.figure(figsize=(8, 8))
    ax = fig.gca()
    contour = ax.contour(array)
    ax.clabel(contour, inline=True, fontsize=10)
    ax.imshow(array[:][:], cmap='coolwarm', alpha=0.7, interpolation='nearest')
    ax.set_title(f'Topology Contour: {refcode} {hkl}')
    plt.show()


def main(args: object) -> None:
    crystal = CrystalReader('CSD').crystal(args.refcode)
    hkl = args.hkl

    surface = particle.Surface(crystal, hkl, surface_size=(2, 2))

    divider = "######" * 5
    output = f"{divider} Surface Analysis {divider}\n"
    output += f"Selected structure          : {surface.identifier}\n"
    output += f"Selected hkl                : {surface.miller_indices}\n"
    output += "\n\n"
    output += f"{divider} Slab Properties {divider}\n"
    output += f"Organic                     : {surface.slab.is_organic}\n"
    output += f"Polymeric                   : {surface.slab.is_polymeric}\n"
    output += f"Number of molecules in slab : {len(surface.slab_molecules)}\n"
    output += f"Slab molecular weight       : {surface.slab.molecular_weight}\n"
    output += f"Number of atoms             : {len(surface.slab.atoms)}\n"
    output += "\n\n"
    output += f"{divider} Topology Properties {divider}\n"
    output += f"Surface area                : {surface.descriptors.surface_area}\n"
    output += f"Rugosity                    : {surface.descriptors.rugosity}\n"
    output += f"RMSD                        : {surface.descriptors.rmsd}\n"
    output += f"Skewness                    : {surface.descriptors.skewness}\n"
    output += f"kurtosis                    : {surface.descriptors.kurtosis}\n"
    output += "\n\n"
    output += f"{divider} Chemical Properties {divider}\n"
    output += f"Density of HB-Donors        : {surface.descriptors.hb_donors}\n"
    output += f"Density of Unsatisfied HB-Donors : {surface.descriptors.unsatisfied_hb_donors}\n"
    output += f"Density of HB-Acceptors     : {surface.descriptors.hb_acceptors}\n"
    output += f"Density of Aromatic Bonds   : {surface.descriptors.aromatic_bonds}\n"
    print(output)
    nodes = list(surface.topology.nodes)

    if args.plot_topology:
        # Rotate nodes to ensure the z-component is perpendicular to surface normal
        normal = surface.miller_indices.plane.normal
        normal = np.array([normal.x, normal.y, normal.z])
        nodes = _rotate_to_xy(normal, nodes)

        x_t = np.array([cor[0] for cor in nodes])
        y_t = np.array([cor[1] for cor in nodes])

        print(f"{divider} Topology Properties {divider}")
        print(f"Dimensions (x,y)              : ({round(x_t.max() - x_t.min())}, {round(y_t.max() - y_t.min())})")
        cell_limits = _node_limits(nodes)
        array = _generate_interpolated_grid_data(nodes, cell_limits,
                                                 grid_spacing=0.2)
        array = gaussian_filter(array, sigma=1.5)
        plot_contour(array, args.refcode, args.hkl)
    else:
        # generates output file
        if args.node_data:
            output += f"{divider} Nodes {divider}\n"
            for node in nodes:
                output += str(node) + "\n"
        with open(os.path.join(os.getcwd(),
                               '{}_{}.txt'.format(args.refcode,
                                                  ''.join(map(str, args.hkl)))),
                  'w') as fh:
            fh.write(output)


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

if __name__ == '__main__':
    args = parse_arguments()
    main(args)