diff --git a/.gitignore b/.gitignore
index b6e4761..bd6ad26 100644
--- a/.gitignore
+++ b/.gitignore
@@ -127,3 +127,6 @@ dmypy.json
# Pyre type checker
.pyre/
+
+# PyCharm
+.idea/
diff --git a/README.md b/README.md
index 70e5607..18725e5 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,57 @@
-# csd-python-api-scripts
-Example scripts using the CSD Python API
+
+
+# CSD Python API Scripts
+
+
+A repository containing scripts that have been created to leverage the toolkit found within
+the [CCDC portfolio](https://www.ccdc.cam.ac.uk/solutions/) that are accessible via
+the [CSD-Python API ](https://www.ccdc.cam.ac.uk/solutions/csd-core/components/csd-python-api/).
+
+The purpose of this platform is to distribute knowledge and allow for scientific collaborations. Scripts are provided on an as-is basis and while their use is not supported we do welcome feedback on potential improvements. All scripts are tested against the latest version of the CSD Python API as installed with the CSD Portfolio.
+
+> For feedback or to report any issues please contact [support@ccdc.cam.ac.uk](support@ccdc.cam.ac.uk)
+
+## Content
+
+```graphql
+.
+├─ assets # Images for documentation
+├─ scripts # Python scripts
+├─ LICENSE
+└─ README.md
+```
+
+## Licensing Requirements
+
+Valid CSD Portfolio licence is required for all scripts, some will require additional licence levels for use, such as Discovery and Materials. Requirements are clearly displayed in the `ReadMe` that accompanies each script.
+
+Don’t have a licence? Contact us [here to request a quote or demonstration.](https://www.ccdc.cam.ac.uk/theccdcprofile/contactus/)
+
+## Downloading Scripts
+
+This section will focus on how to download multiple or individual scripts:
+
+### 1. Downloading All Scripts
+
+
+#### Downloading Zip file:
+Click `Code` and Select `Download Zip`:
+
+
+
+
+### 2. Downloading Individual Scripts
+
+#### TODO - Add gif of downloading single file
+#### TODO - write instructions on how to download
+
+
+## Running scripts through the CSD Python API Miniconda installed
+
+This requires you to already have the CSD Python API installed.
+
+Windows:
+```cmd
+"\CCDC\Python_API_2022\miniconda\Scripts\activate.bat"
+python script_example.py
+```
\ No newline at end of file
diff --git a/assets/csd-python-api-logo.png b/assets/csd-python-api-logo.png
new file mode 100644
index 0000000..c4f1242
Binary files /dev/null and b/assets/csd-python-api-logo.png differ
diff --git a/assets/download_zip.png b/assets/download_zip.png
new file mode 100644
index 0000000..3250716
Binary files /dev/null and b/assets/download_zip.png differ
diff --git a/scripts/ReadMe.md b/scripts/ReadMe.md
new file mode 100644
index 0000000..a98439e
--- /dev/null
+++ b/scripts/ReadMe.md
@@ -0,0 +1,16 @@
+## Contents
+
+This folder contains scripts submitted by users or CCDC scientists for anyone to use freely.
+
+### Hydrogen bond propensity :
+- Writes a .docx report of a hydrogen bond propensity calculation for any given mol2/refcode.
+
+### Multi-component hydrogen bond propensity:
+- Performs a multi-component HBP calculation for a given library of co-formers.
+
+### Packing similarity dendrogram:
+- Construct a dendrogram for an input set of structures based on packing-similarity
+analysis
+
+## Tips
+### Searching tips:
\ No newline at end of file
diff --git a/scripts/hydrogen_bond_propensity/ReadMe.md b/scripts/hydrogen_bond_propensity/ReadMe.md
new file mode 100644
index 0000000..58f07b9
--- /dev/null
+++ b/scripts/hydrogen_bond_propensity/ReadMe.md
@@ -0,0 +1,51 @@
+# Hydrogen Bond Propensity
+
+Writes a .docx report of a hydrogen bond propensity calculation
+
+----
+Included in the report:
+- Hydrogen bond propensity chart.
+- Table of Predicted intermolecular hydrogen bond propensities.
+- Table of Hydrogen bond coordination likelihood.
+- Table of Hydrogen Bond Donor/Acceptor Definitions
+- Functional Group Definitions for CSD Substructure Searches
+ - Search Results
+ - Number of hits selected per functional group
+ - Hydrogen Bond Analysis of Training Dataset
+ - Hydrogen bond donor and acceptor outcomes
+
+## Example
+
+
+
+
+
+
+
+## Requirements
+
+- ```matplotlib, docxtpl```
+
+- hydrogen_bond_propensity_report.docx
+## Licensing Requirements
+- CSD-Core
+
+## Instructions on running
+
+Using mol2 file
+
+```cmd
+"C:\Program Files\CCDC\Python_API_2022\miniconda\Scripts\activate.bat"
+python hydrogen_bond_propensity_report.py structure.mol2
+```
+
+Options:
+```
+- d | --directory - sets the working directory
+- n | --noopen - Do not automatically open the generated output file.
+```
+## Author
+
+_Andrew Maloney_ - CCDC (2017)
+
+> For feedback or to report any issues please contact [support@ccdc.cam.ac.uk](support@ccdc.cam.ac.uk)
\ No newline at end of file
diff --git a/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_1.png b/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_1.png
new file mode 100644
index 0000000..3724ac0
Binary files /dev/null and b/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_1.png differ
diff --git a/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_2.png b/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_2.png
new file mode 100644
index 0000000..b10f6fc
Binary files /dev/null and b/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_2.png differ
diff --git a/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_3.png b/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_3.png
new file mode 100644
index 0000000..17b998f
Binary files /dev/null and b/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_3.png differ
diff --git a/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_4.png b/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_4.png
new file mode 100644
index 0000000..411a958
Binary files /dev/null and b/scripts/hydrogen_bond_propensity/assets/HXACAN_Report_Screenshot_4.png differ
diff --git a/scripts/hydrogen_bond_propensity/hydrogen_bond_propensity_report.docx b/scripts/hydrogen_bond_propensity/hydrogen_bond_propensity_report.docx
new file mode 100644
index 0000000..9fe5da1
Binary files /dev/null and b/scripts/hydrogen_bond_propensity/hydrogen_bond_propensity_report.docx differ
diff --git a/scripts/hydrogen_bond_propensity/hydrogen_bond_propensity_report.py b/scripts/hydrogen_bond_propensity/hydrogen_bond_propensity_report.py
new file mode 100644
index 0000000..256126e
--- /dev/null
+++ b/scripts/hydrogen_bond_propensity/hydrogen_bond_propensity_report.py
@@ -0,0 +1,348 @@
+#!/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-08-10: Created by Andy Maloney, the Cambridge Crystallographic Data Centre
+# 2020-08-21: made available by the Cambridge Crystallographic Data Centre
+#
+
+"""
+hydrogen_bond_propensity_report.py
+- Writes a .docx report of a hydrogen bond propensity calculation
+"""
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+
+from ccdc import io
+from ccdc.diagram import DiagramGenerator
+from ccdc.search import SubstructureSearch, ConnserSubstructure
+from ccdc.descriptors import CrystalDescriptors
+
+import argparse
+import os
+import sys
+import subprocess
+import csv
+
+try:
+ import warnings
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore", category=DeprecationWarning)
+ import docxtpl
+ from docx.shared import Cm
+except ImportError:
+ error_message = """
+ The python-docx-template templating engine needed by this script could not
+ be found. Please run "{} -m pip install docxtpl"
+ to try to fix the issue.\nYou may need administrator's rights to do this.
+ """.format(sys.executable)
+ raise ImportError(error_message)
+
+SCRIPT_DIR = os.path.dirname(__file__)
+TEMPLATE_FILENAME = 'hydrogen_bond_propensity_report.docx'
+TEMPLATE_FILE = os.path.join(SCRIPT_DIR, TEMPLATE_FILENAME)
+
+
+###############################################################################
+def make_diagram(mol, directory):
+ # Generates a diagram from a given structure
+ molecule_diagram_generator = DiagramGenerator()
+ molecule_diagram_generator.settings.line_width = 1.6
+ molecule_diagram_generator.settings.font_size = 12
+ molecule_diagram_generator.settings.image_height = 300
+ img = molecule_diagram_generator.image(mol)
+ fname = str(os.path.join(directory, '%s_diagram.png' % mol.identifier))
+ if img:
+ img.save(fname)
+ return fname
+
+
+def fg_diagram(mol, directory, con):
+ # Create highlighted functional group diagrams
+ diagram_generator = DiagramGenerator()
+ diagram_generator.settings.shrink_symbols = False
+ diagram_generator.settings.element_coloring = False
+ searcher = SubstructureSearch()
+ searcher.add_substructure(ConnserSubstructure(os.path.join(directory, "%s.con" % con)))
+ hits = searcher.search(mol)
+ selection = hits[0].match_atoms()
+ img = diagram_generator.image(mol, highlight_atoms=selection)
+ fname = str(os.path.join(directory, '%s.png' % con))
+ if img:
+ img.save(fname)
+ return fname
+
+
+def add_picture_subdoc(picture_location, docx_template, cm=7):
+ # This function adds a picture to the .docx file
+ return docxtpl.InlineImage(
+ docx_template, image_descriptor=picture_location, width=Cm(cm))
+
+
+def launch_word_processor(output_file):
+ """Open the default application for output_file across platforms."""
+ if sys.platform == 'win32':
+ os.startfile(output_file)
+ elif sys.platform.startswith('linux'):
+ subprocess.Popen(['xdg-open', output_file])
+ else:
+ subprocess.Popen(['open', output_file])
+
+
+def propensity_calc(crystal, directory):
+ # Perform a Hydrogen Bond Propensity calculation
+
+ # Provide settings for the calculation
+ settings = CrystalDescriptors.HBondPropensities.Settings()
+ settings.working_directory = directory
+ settings.hbond_criterion.require_hydrogens = True
+ settings.hbond_criterion.path_length_range = (3, 999)
+
+ # Set up the HBP calculator
+ hbp = CrystalDescriptors.HBondPropensities(settings)
+
+ # Set up the target structure for the calculation
+ hbp.set_target(crystal)
+
+ print(hbp.functional_groups)
+
+ # Generate Training Dataset
+
+ hbp.match_fitting_data(count=300) # set to >300
+
+ hbp.analyse_fitting_data()
+
+ for d in hbp.donors:
+ print(d.identifier, d.npositive, d.nnegative)
+ for a in hbp.acceptors:
+ print(a.identifier, a.npositive, a.nnegative)
+
+ # Perform regression
+ model = hbp.perform_regression()
+
+ print(model.equation)
+ print('Area under ROC curve: {} -- {}'.format(round(model.area_under_roc_curve, 3), model.advice_comment))
+
+ propensities = hbp.calculate_propensities()
+ if len(hbp.intra_propensities) > 0:
+ intra_flag = True
+ else:
+ intra_flag = False
+ groups = hbp.generate_hbond_groupings()
+ observed_group = hbp.target_hbond_grouping()
+
+ return hbp.functional_groups, hbp.fitting_data, hbp.donors, hbp.acceptors, model, \
+ propensities, intra_flag, groups, observed_group
+
+
+def coordination_scores_calc(crystal, directory):
+ # Calculate coordination scores for the target structure
+
+ # Provide settings for the calculation
+ settings = CrystalDescriptors.HBondCoordination.Settings()
+ settings.working_directory = directory
+
+ # Set up the coordination scores calculator
+ coordination_calc = CrystalDescriptors.HBondCoordination(settings)
+ predictions = coordination_calc.predict(crystal)
+ return predictions
+
+
+def format_scores(scores, das, d_type):
+ # Reformat the coordination scores to make report writing easier
+ formatted_scores = {}
+ for da in das:
+ preds = scores.predictions_for_label(da.label, d_type)[1]
+ formatted_scores[da.label] = preds
+ return formatted_scores
+
+
+def normalize_molecule(molecule):
+ # Normalise bond types for the input structure (important for cifs)
+ molecule.assign_bond_types(which='unknown')
+ molecule.standardise_aromatic_bonds()
+ molecule.standardise_delocalised_bonds()
+
+
+def chart_output(groups, work_directory, structure):
+ # Write out the data points of the HBP chart to a file
+ with open(os.path.join(work_directory, '%s_chart_data.csv' % structure), 'w') as outfile:
+ csv_writer = csv.writer(outfile)
+ csv_writer.writerow(['Mean Propensity', 'Mean Coordination Score', 'Hydrogen Bonds'])
+ for group in groups:
+ csv_writer.writerow(
+ [group.hbond_score,
+ group.coordination_score,
+ '; '.join(['%s - %s' % (g.donor.label, g.acceptor.label) for g in group.hbonds])]
+ )
+
+
+def main(structure, directory, csdrefcode, noopen=False):
+ # This looks for the .docx template that is used to generate the report from
+ if os.path.isfile(TEMPLATE_FILE):
+ docx_template = docxtpl.DocxTemplate(TEMPLATE_FILE)
+ else:
+ print('Error! {} not found!'.format(TEMPLATE_FILENAME))
+ quit()
+
+ # This loads up the CSD if a refcode is requested, otherwise loads the structural file supplied
+ if csdrefcode:
+ try:
+ crystal = io.CrystalReader('CSD').crystal(str(structure))
+ except RuntimeError:
+ print('Error! %s is not in the database!' % str(structure))
+ quit()
+ else:
+ crystal = io.CrystalReader(str(structure))[0]
+
+ # If there are atoms without sites, or there are no atoms in the structure, then HBP cannot be performed
+ molecule = crystal.molecule
+ if not molecule.all_atoms_have_sites or len(molecule.atoms) == 0:
+ print('Error! Not all atoms in %s have sites!' % str(structure))
+ quit()
+
+ # Bond types need to be standardised
+ normalize_molecule(molecule)
+ crystal.molecule = molecule
+
+ # Set up a work directory for the HBP files
+ work_directory = os.path.join(directory, str(structure).split('.')[0])
+
+ # Get all the necessary data from a HBP calculation
+ functional_groups, fitting_data, donors, acceptors, model, propensities, intra_flag,\
+ groups, observed_groups = propensity_calc(crystal, work_directory)
+
+ # Calculate the coordination scores separately
+ coordination_scores = coordination_scores_calc(crystal, work_directory)
+
+ # Create the HBP chart output as a separate file
+ chart_output(groups, work_directory, structure)
+
+ # Set up some dictionaries and fill them with coordination score data, as well as extra information for the report
+ dscores = {}
+ ascores = {}
+ dcoord_bg = {}
+ acoord_bg = {}
+
+ for d in donors:
+ coord_cols = [i for i in range(len(coordination_scores.predictions_for_label(d.label, 'd')[1]))]
+ dscores[d.label] = [round((coordination_scores.predictions_for_label(d.label, 'd')[1])[j], 3)
+ for j in coord_cols]
+ dcoord_bg[d.label] = ['FFFFFF' for j in coord_cols]
+ coord_value = coordination_scores.predictions_for_label(d.label, 'd')[0]
+ if coord_value == dscores[d.label].index(max(dscores[d.label])):
+ dcoord_bg[d.label][coord_value] = '7FFF00'
+ else:
+ dcoord_bg[d.label][coord_value] = 'FF0000'
+
+ for a in acceptors:
+ ascores[a.label] = [round((coordination_scores.predictions_for_label(a.label, 'a')[1])[k], 3)
+ for k in coord_cols]
+ acoord_bg[a.label] = ['FFFFFF' for j in coord_cols]
+ coord_value = coordination_scores.predictions_for_label(a.label, 'a')[0]
+ if coord_value == ascores[a.label].index(max(ascores[a.label])):
+ acoord_bg[a.label][coord_value] = '7FFF00'
+ else:
+ acoord_bg[a.label][coord_value] = 'FF0000'
+
+ # Generate more information for the report
+ coefficients = model.coefficients
+ don = list(set(list("%s_d" % p.donor_label.split(" ")[0] for p in propensities)))
+ acc = list(set(list("%s_a" % p.acceptor_label.split(" ")[0] for p in propensities)))
+
+ # Generate the HBP chart
+ figure = plt.scatter([group.hbond_score for group in groups],
+ [group.coordination_score for group in groups],
+ s=50)
+ ax2 = plt.scatter(observed_groups.hbond_score,
+ observed_groups.coordination_score,
+ c='red', marker='*', s=250)
+ plt.origin = 'upper'
+ plt.xlim(0, 1.0)
+ plt.ylim(-1.0, 0)
+ 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')
+ figure_location = os.path.join(directory, '%s_standard_chart.png' % crystal.identifier)
+ plt.savefig(figure_location)
+
+ chart = add_picture_subdoc(figure_location, docx_template, cm=16)
+
+ diagram_file = make_diagram(crystal.molecule, directory)
+ diagram = add_picture_subdoc(diagram_file, docx_template)
+
+ con_files = [f[:-4] for f in os.listdir(work_directory) if f.endswith(".con")]
+ fg_pics = [fg_diagram(crystal.molecule, work_directory, con) for con in con_files]
+ fg_diagrams = {con: add_picture_subdoc(os.path.join(work_directory, '%s.png' % con), docx_template)
+ for con in con_files}
+
+ # The context is the information that is given to the template to allow it to be populated
+ context = {
+ # Title page
+ 'identifier': crystal.identifier,
+ 'chart': chart,
+ 'propensities': propensities,
+ 'intra_flag': intra_flag,
+ 'coord_cols': coord_cols,
+ 'donors': donors,
+ 'acceptors': acceptors,
+ 'dscores': dscores,
+ 'dbg': dcoord_bg,
+ 'ascores': ascores,
+ 'abg': acoord_bg,
+ 'diagram': diagram,
+ 'don': don,
+ 'acc': acc,
+ 'fg_diagrams': fg_diagrams,
+ 'functional_groups': functional_groups,
+ 'data': fitting_data,
+ 'len_data': len(fitting_data),
+ 'coefficients': coefficients,
+ 'model': model,
+ }
+
+ # Send all the information to the template file then open up the final report
+ docx_template.render(context)
+ output_file = os.path.join(directory, '%s_propensity_report.docx' % crystal.identifier)
+ docx_template.save(output_file)
+
+ if not noopen:
+ launch_word_processor(output_file)
+ print('Output file written to %s' % output_file)
+
+
+if __name__ == '__main__':
+ # Set up the necessary arguments to run the script
+ parser = argparse.ArgumentParser(
+ formatter_class=argparse.RawDescriptionHelpFormatter,
+ description=__doc__)
+ parser.add_argument('input_structure', type=str,
+ help='Refcode or mol2 file of the component to be screened')
+ parser.add_argument('-d', '--directory', default=os.getcwd(),
+ help='the working directory for the calculation')
+ parser.add_argument('-n', '--noopen', action='store_true',
+ help='Do not automatically open the generated output file.')
+
+ args = parser.parse_args()
+
+ refcode = False
+
+ if not os.path.isfile(args.input_structure):
+ if len(str(args.input_structure).split('.')) == 1:
+ refcode = True
+ else:
+ parser.error('%s - file not found.' % args.input_structure)
+ if not refcode:
+ args.directory = os.path.dirname(os.path.abspath(args.input_structure))
+ elif not os.path.isdir(args.directory):
+ os.makedirs(args.directory)
+
+ main(args.input_structure, args.directory, refcode, args.noopen)
diff --git a/scripts/multi_component_hydrogen_bond_propensity/ReadMe.md b/scripts/multi_component_hydrogen_bond_propensity/ReadMe.md
new file mode 100644
index 0000000..39aed38
--- /dev/null
+++ b/scripts/multi_component_hydrogen_bond_propensity/ReadMe.md
@@ -0,0 +1,51 @@
+# Multi-Component Hydrogen Bond Propensity
+
+## Summary
+
+Performs a multi-component Hydrogen bond propensity calculation for a given library of co-formers
+## Example
+
+Individual reports are generated for each coformer stored in indvidual folders. A summary multicomponent report with rankings is also generated.
+
+Individual reports include:
+
+- Predicted intermolecular hydrogen bond propensities
+- Hydrogen bond coordination likelihood
+
+Summary report includes:
+- Chemical Diagram
+- Table of ranked components
+- Multi-component hydrogen-bond propensity chart
+
+## Requirements
+- docxtpl
+- multi_component_hydrogen_bond_propensity_report.docx
+- multi_component_pair_hbp_report.docx
+
+## Licensing Requirements
+
+- CSD-Core
+
+## Instructions on running
+
+```cmd
+python multi_component_hydrogen_bond_propensity_report.py HXACAN28
+```
+
+```cmd
+positional arguments:
+ input_structure Refcode or mol2 file of the component to be screened
+
+optional arguments:
+ -h, --help show this help message and exit
+ -d DIRECTORY, --directory DIRECTORY
+ the working directory for the calculation
+ -c COFORMER_LIBRARY, --coformer_library COFORMER_LIBRARY
+ the directory of the desired coformer library
+```
+
+The default coformer library is the one supplied with your Mercury install in ```\Mercury\molecular_libraries\ccdc_coformers```
+## Author
+_Andrew Maloney_ (CCDC) 2017
+
+> For feedback or to report any issues please contact [support@ccdc.cam.ac.uk](support@ccdc.cam.ac.uk)
\ No newline at end of file
diff --git a/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.docx b/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.docx
new file mode 100644
index 0000000..077b264
Binary files /dev/null and b/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.docx differ
diff --git a/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.py b/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.py
new file mode 100644
index 0000000..8486e2e
--- /dev/null
+++ b/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.py
@@ -0,0 +1,404 @@
+#!/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-08-10: Created by Andy Maloney, the Cambridge Crystallographic Data Centre
+# 2020-08-21: made available by the Cambridge Crystallographic Data Centre
+#
+"""
+multi_component_hydrogen_bond_propensity_report.py
+ - Performs a multi-component HBP calculation for a given library of co-formers
+"""
+
+import matplotlib
+matplotlib.use('TkAgg')
+import matplotlib.pyplot as plt
+
+from ccdc import io, molecule
+from ccdc.diagram import DiagramGenerator
+from ccdc.descriptors import CrystalDescriptors
+
+import sys
+import os
+import glob
+import argparse
+import tempfile
+import subprocess
+
+try:
+ import warnings
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore", category=DeprecationWarning)
+ import docxtpl
+ from docx.shared import Cm
+except ImportError:
+ error_message = """
+ The python-docx-template templating engine needed by this script could not
+ be found. Please run "{} -m pip install docxtpl"
+ to try to fix the issue.\nYou may need administrator's rights to do this.
+ """.format(sys.executable)
+ raise ImportError(error_message)
+
+SCRIPT_DIR = os.path.dirname(__file__)
+TEMPLATE_FILENAME = 'multi_component_hydrogen_bond_propensity_report.docx'
+TEMPLATE_FILE = os.path.join(SCRIPT_DIR, TEMPLATE_FILENAME)
+PAIR_TEMPLATE_FILENAME = 'multi_component_pair_hbp_report.docx'
+PAIR_TEMPLATE_FILE = os.path.join(SCRIPT_DIR, PAIR_TEMPLATE_FILENAME)
+
+
+###############################################################################
+def cm2inch(*tupl):
+ inch = 2.54
+ if isinstance(tupl[0], tuple):
+ return tuple(i/inch for i in tupl[0])
+ else:
+ return tuple(i/inch for i in tupl)
+
+
+def launch_word_processor(output_file):
+ '''This function launches the platform specific word processor
+ when not running under continuous integration'''
+ if 'TEAMCITY_VERSION' in os.environ:
+ return
+ if sys.platform == 'win32':
+ os.startfile(output_file)
+ elif sys.platform.startswith('linux'):
+ subprocess.Popen(['xdg-open', output_file])
+ else:
+ subprocess.Popen(['open', output_file])
+
+
+def make_diagram(mol, directory):
+ # Generates a diagram from a given structure
+ molecule_diagram_generator = DiagramGenerator()
+ molecule_diagram_generator.settings.line_width = 1.6
+ molecule_diagram_generator.settings.font_size = 12
+ molecule_diagram_generator.settings.image_height = 300
+ img = molecule_diagram_generator.image(mol)
+ fname = str(os.path.join(directory, '%s_diagram.png' % mol.identifier))
+ if img:
+ img.save(fname)
+ return fname
+
+
+def make_mc_chart(dictionary, directory, mol):
+ results = [value[0] for key, value in dictionary if isinstance(value[0], float)]
+ ymin = min(results)
+ ymax = max(results)
+ indices = range(1, len(results) + 1)
+
+ fig = plt.figure(figsize=cm2inch(22, 18))
+ ax = fig.add_subplot(1, 1, 1)
+ color = 'cornflowerblue'
+ color1 = 'royalblue'
+ ls = ''
+ plt.plot(indices, results, marker='D', markersize=10, color=color, ls=ls, markeredgecolor=color1, alpha=0.7)
+ plt.axhline(y=0, color='gray')
+ plt.xlabel('Co-Former Rank', fontweight='bold', fontsize='12')
+ plt.ylabel('Multi-Component Score', fontweight='bold', fontsize='12')
+ plt.title('MCHBP screening results', fontweight='bold', fontsize='12')
+ ax.axhspan(0.025, ymax + 0.1, facecolor='lightgreen', alpha=0.5)
+ ax.axhspan(-0.025, ymin - 0.1, facecolor='lightpink', alpha=0.5)
+ ax.axhspan(0.025, -0.025, facecolor='grey', alpha=0.5)
+ plt.ylim(ymin - 0.025, ymax + 0.025)
+ fname = str(os.path.join(directory, '%s_MC_HBP_plot.png' % mol.identifier))
+ plt.savefig(fname, format='png', dpi=600)
+ return fname
+
+
+def add_picture_subdoc(picture_location, docx_template, wd=7):
+ # This function adds a picture to the .docx file
+ return docxtpl.InlineImage(
+ docx_template, image_descriptor=picture_location, width=Cm(wd))
+
+
+def propensity_calc(crystal, directory):
+ # Perform a Hydrogen Bond Propensity calculation
+
+ # Provide settings for the calculation
+ settings = CrystalDescriptors.HBondPropensities.Settings()
+ settings.working_directory = directory
+ settings.hbond_criterion.require_hydrogens = True
+ settings.hbond_criterion.path_length_range = (3, 999)
+
+ # Set up the HBP calculator
+ hbp = CrystalDescriptors.HBondPropensities(settings)
+
+ # Set up the target structure for the calculation
+ hbp.set_target(crystal)
+
+ print(hbp.functional_groups)
+
+ # Generate Training Dataset
+
+ hbp.match_fitting_data(count=300) # set to >300
+
+ hbp.analyse_fitting_data()
+
+ for d in hbp.donors:
+ print(d.identifier, d.npositive, d.nnegative)
+ for a in hbp.acceptors:
+ print(a.identifier, a.npositive, a.nnegative)
+
+ # Perform regression
+ model = hbp.perform_regression()
+ print(model.equation)
+ print('Area under ROC curve: {} -- {}'.format(round(model.area_under_roc_curve, 3), model.advice_comment))
+ hbp.calculate_propensities()
+ propensities = hbp.inter_propensities
+
+ return propensities, hbp.donors, hbp.acceptors
+
+
+def coordination_scores_calc(crystal, directory):
+ # Calculate coordination scores for the target structure
+
+ # Provide settings for the calculation
+ settings = CrystalDescriptors.HBondCoordination.Settings()
+ settings.working_directory = directory
+ settings.hbond_criterion.require_hydrogens = True
+ settings.hbond_criterion.path_length_range = (3, 999)
+
+ # Set up the coordination scores calculator
+ coordination_calc = CrystalDescriptors.HBondCoordination(settings)
+ predictions = coordination_calc.predict(crystal)
+ return predictions
+
+
+def format_scores(scores, das, d_type):
+ # Reformat the coordination scores to make report writing easier
+ formatted_scores = {}
+ for da in das:
+ preds = scores.predictions_for_label(da.label, d_type)[1]
+ formatted_scores[da.label] = preds
+ return formatted_scores
+
+
+def get_mc_scores(propensities, identifier):
+ # Calculates the multi-component scores from the individual HBP calculation
+ AA_propensities = []
+ BB_propensities = []
+ AB_propensities = []
+ BA_propensities = []
+
+ for p in propensities:
+ t = "%s_d" % p.donor_label.split(" ")[0], "%s_a" % p.acceptor_label.split(" ")[0]
+ if '_A_' in t[0] and '_A_' in t[1]:
+ AA_propensities.append(p.propensity)
+ elif '_B_' in t[0] and '_B_' in t[1]:
+ BB_propensities.append(p.propensity)
+ elif '_A_' in t[0] and '_B_' in t[1]:
+ AB_propensities.append(p.propensity)
+ elif '_B_' in t[0] and '_A_' in t[1]:
+ BA_propensities.append(p.propensity)
+ max_AA = max(AA_propensities) if len(AA_propensities) > 0 else 0.0
+ max_BB = max(BB_propensities) if len(BB_propensities) > 0 else 0.0
+ max_AB = max(AB_propensities) if len(AB_propensities) > 0 else 0.0
+ max_BA = max(BA_propensities) if len(BA_propensities) > 0 else 0.0
+ max_list = [max_AA, max_BB, max_AB, max_BA]
+ max_keys = ['A:A', 'B:B', 'A:B', 'B:A']
+ max_mc = max(max_list[2], max_list[3])
+ max_sc = max(max_list[0], max_list[1])
+
+ return [round((max_mc - max_sc), 2),
+ max_keys[max_list.index(max(max_list))],
+ round(max_mc, 2),
+ round(max_list[0], 2),
+ round(max_list[1], 2),
+ identifier]
+
+
+def make_pair_file(api_molecule, tempdir, f, i):
+ # Creates a file for the api/coformer pair
+ with io.MoleculeReader(f) as reader:
+ coformer_molecule = reader[0]
+ coformer_name = coformer_molecule.identifier
+ molecule_pair = make_molecule_pair(api_molecule, coformer_molecule, i)
+ molecule_file = os.path.join(tempdir, '%s.mol2' % molecule_pair.identifier)
+ with io.MoleculeWriter(molecule_file) as writer:
+ writer.write(molecule_pair)
+ return molecule_file, coformer_name
+
+
+def make_molecule_pair(api_molecule, coformer_molecule, i):
+ # Creates the multi-component system for each api/coformer pair
+ new_file_name = '%s_%d' % (api_molecule.identifier, i)
+ molecule_pair = molecule.Molecule(new_file_name)
+ molecule_pair.add_molecule(api_molecule)
+ molecule_pair.add_molecule(coformer_molecule)
+ molecule_pair.normalise_labels()
+
+ atom_list = []
+ for atom in molecule_pair.components[0].atoms:
+ atom_list.append(atom.label)
+
+ for atom in molecule_pair.atoms:
+ if atom.label in atom_list:
+ atom.label += '_A'
+ else:
+ atom.label += '_B'
+
+ return molecule_pair
+
+
+def pair_output(identifier, propensities, donors, acceptors, coordination_scores, directory):
+ # Writes out the output from a single HBP calculation for the multi-component pair
+
+ # This looks for the .docx template that is used to generate the report from
+ if os.path.isfile(PAIR_TEMPLATE_FILE):
+ docx_template = docxtpl.DocxTemplate(PAIR_TEMPLATE_FILE)
+ else:
+ print('Error! {} not found!'.format(PAIR_TEMPLATE_FILENAME))
+ quit()
+
+ dscores = {}
+ ascores = {}
+
+ for d in donors:
+ coord_cols = [i for i in range(len(coordination_scores.predictions_for_label(d.label, 'd')[1]))]
+ dscores[d.label] = [round((coordination_scores.predictions_for_label(d.label, 'd')[1])[j], 3)
+ for j in coord_cols]
+
+ for a in acceptors:
+ ascores[a.label] = [round((coordination_scores.predictions_for_label(a.label, 'a')[1])[k], 3)
+ for k in coord_cols]
+
+ context = {
+ 'identifier': identifier,
+ 'propensities': propensities,
+ 'coord_cols': coord_cols,
+ 'donors': donors,
+ 'dscores': dscores,
+ 'acceptors': acceptors,
+ 'ascores': ascores
+ }
+ docx_template.render(context)
+ output_file = os.path.join(directory, '%s_pair_output.docx' % identifier)
+ docx_template.save(output_file)
+
+
+def make_mc_report(identifier, results, directory, diagram_file, chart_file):
+ # Write the MC-HBP report from the results
+
+ # This looks for the .docx template that is used to generate the report from
+ if os.path.isfile(TEMPLATE_FILE):
+ docx_template = docxtpl.DocxTemplate(TEMPLATE_FILE)
+ else:
+ print('Error! {} not found!'.format(TEMPLATE_FILENAME))
+ quit()
+
+ # Generate content for the report
+ diagram = add_picture_subdoc(diagram_file, docx_template)
+ chart = add_picture_subdoc(chart_file, docx_template,wd=18)
+
+ # The context is the information that is given to the template to allow it to be populated
+ context = {
+ 'identifier': str(identifier).split('.')[0],
+ 'diagram': diagram,
+ 'chart': chart,
+ 'results': results
+ }
+
+ # Send all the information to the template file then open up the final report
+ docx_template.render(context)
+ output_file = os.path.join(directory, '%s_MC_HBP_report.docx' % str(identifier).split('.')[0])
+ docx_template.save(output_file)
+
+ launch_word_processor(output_file)
+
+
+def main(structure, work_directory, library, csdrefcode):
+ # This loads up the CSD if a refcode is requested, otherwise loads the structural file supplied
+ if csdrefcode:
+ try:
+ crystal = io.CrystalReader('CSD').crystal(structure)
+ except RuntimeError:
+ print('Error! %s is not in the database!' % structure)
+ quit()
+ else:
+ crystal = io.CrystalReader(structure)[0]
+
+ # If there are atoms without sites, or there are no atoms in the structure, then HBP cannot be performed
+ api_molecule = crystal.molecule
+ if not api_molecule.all_atoms_have_sites or len(api_molecule.atoms) == 0:
+ print('Error! Not all atoms in %s have sites!' % structure)
+ quit()
+
+ # find the coformers and set up the calculations
+ coformer_files = glob.glob(os.path.join(library, '*.mol2'))
+ tempdir = tempfile.mkdtemp()
+ mc_dictionary = {}
+
+ # for each coformer in the library, make a pair file for the api/coformer and run a HBP calculation
+ for i, f in enumerate(coformer_files):
+ molecule_file, coformer_name = make_pair_file(api_molecule, tempdir, f, i + 1)
+ print(coformer_name)
+ crystal_reader = io.CrystalReader(molecule_file)
+ crystal = crystal_reader[0]
+
+ directory = os.path.join(os.path.abspath(work_directory), crystal.identifier)
+
+ try:
+ propensities, donors, acceptors = propensity_calc(crystal, directory)
+ coordination_scores = coordination_scores_calc(crystal, directory)
+ pair_output(crystal.identifier, propensities, donors, acceptors, coordination_scores, directory)
+ mc_dictionary[coformer_name] = get_mc_scores(propensities, crystal.identifier)
+
+ except RuntimeError:
+ print("Propensity calculation failure for %s!" % coformer_name)
+ mc_dictionary[coformer_name] = ["N/A", "N/A", "N/A", "N/A", "N/A", crystal.identifier]
+
+ # Make sense of the outputs of all the calculations
+ mc_hbp_screen = sorted(mc_dictionary.items(), key=lambda e: e[1][0], reverse=True)
+ diagram_file = make_diagram(api_molecule, work_directory)
+ chart_file = make_mc_chart(mc_hbp_screen, directory, api_molecule)
+ make_mc_report(structure, mc_hbp_screen, work_directory, diagram_file, chart_file)
+
+
+if __name__ == '__main__':
+ # Set up the necessary arguments to run the script
+ if sys.platform == 'win32':
+ ccdc_coformers_dir = os.path.join(
+ os.path.dirname(io.csd_directory()),
+ 'Mercury',
+ 'molecular_libraries',
+ 'ccdc_coformers'
+ )
+ else:
+ ccdc_coformers_dir = os.path.join(
+ os.path.dirname(io.csd_directory()),
+ 'molecular_libraries',
+ 'ccdc_coformers'
+ )
+ parser = argparse.ArgumentParser(
+ formatter_class=argparse.RawDescriptionHelpFormatter,
+ description=__doc__)
+ parser.add_argument('input_structure', type=str,
+ help='Refcode or mol2 file of the component to be screened')
+ parser.add_argument('-d', '--directory', default=os.getcwd(),
+ help='the working directory for the calculation')
+ parser.add_argument('-c', '--coformer_library', type=str,
+ help='the directory of the desired coformer library',
+ default=ccdc_coformers_dir)
+
+ args = parser.parse_args()
+
+ refcode = False
+
+ if not os.path.isfile(args.input_structure):
+ if len(str(args.input_structure).split('.')) == 1:
+ refcode = True
+ else:
+ parser.error('%s - file not found.' % args.input_structure)
+ if not refcode:
+ args.directory = os.path.dirname(os.path.abspath(args.input_structure))
+ elif not os.path.isdir(args.directory):
+ os.makedirs(args.directory)
+ if not os.path.isdir(args.coformer_library):
+ parser.error('%s - library not found.' % args.coformer_library)
+
+ main(args.input_structure, args.directory, args.coformer_library, refcode)
diff --git a/scripts/multi_component_hydrogen_bond_propensity/multi_component_pair_hbp_report.docx b/scripts/multi_component_hydrogen_bond_propensity/multi_component_pair_hbp_report.docx
new file mode 100644
index 0000000..5de0323
Binary files /dev/null and b/scripts/multi_component_hydrogen_bond_propensity/multi_component_pair_hbp_report.docx differ
diff --git a/scripts/new_script_readme_template/ReadMe.md b/scripts/new_script_readme_template/ReadMe.md
new file mode 100644
index 0000000..963d76e
--- /dev/null
+++ b/scripts/new_script_readme_template/ReadMe.md
@@ -0,0 +1,37 @@
+# Script Name
+
+----
+
+## Summary
+
+\#Description of what the script does
+
+\#Example output showing what the user can expect to see
+
+## Requirements
+
+\#Minimum CSD Python API version
+
+\#List of additional required CCDC modules, e.g. ccdc_rp
+
+\#List of additional non-CCDC python modules required, e.g. docxtpl
+
+\#Any additional files/resources required, e.g. additional files from this repository script folder
+
+## Licensing Requirements
+
+\#CCDC Licence requirements, e.g. CSD-Core, CSD-Materials, CSD-Enterprise, RP
+
+## Instructions on running
+
+\#Command line usage
+
+\#Arguments available
+
+\#Additional resources required, e.g. file from the script repository and how that is used (e.g. in same folder as python script, used in an argument etc.)
+
+## Author
+
+\#Author name and date created
+
+> For feedback or to report any issues please contact [support@ccdc.cam.ac.uk](support@ccdc.cam.ac.uk)
\ No newline at end of file
diff --git a/scripts/packing_similarity_dendrogram/ReadMe.md b/scripts/packing_similarity_dendrogram/ReadMe.md
new file mode 100644
index 0000000..40b9384
--- /dev/null
+++ b/scripts/packing_similarity_dendrogram/ReadMe.md
@@ -0,0 +1,167 @@
+# Packing_Similarity_Dendrogram.py
+
+## Summary:
+
+This command-line script is designed to compare the packing of a set of input structures of a molecule (polymorphs,
+co-crystals, solvates, and hydrates). An all-to-all comparison of the structures is performed, considering _only the
+heaviest components_ in each structure, and a packing-similarity dendrogram or tree is constructed using hierarchical
+clustering. This dendrogram shows the similarity between groups of structures and how these groups relate to one
+another.
+
+## Requirements:
+
+CSD Python API (v. 1.0 or later), `matplotlib` and standard python packages.
+
+> Instructions on running example scripts through the CSD Python API miniconda installation can be found in [the repository ReadMe.md](https://github.com/ccdc-opensource/csd-python-api-scripts#readme)
+
+## Licensing Requirements
+
+- CSD-Materials or CSD-Enterprise
+
+## Usage and Output:
+
+`python Packing_Similarity_Dendrogram.py -h` will show the help text:
+
+```
+usage: Packing_Similarity_Dendrogram.py [-h] [-m similarity_matrix.txt]
+[-ns 25] [-nm 15] [-o]
+[--allow_molecular_differences]
+[--clustering_type {complete,single,average}]
+[-s] [-ct 0.5] [-at 25] [-dt 0.25]
+input_file
+
+Packing_Similarity_Dendrogram.py - Construct a dendrogram for an input set of structures based
+on packing-similarity analysis
+
+positional arguments:
+input_file Set of structures to perform analysis on
+[.mol2/cif/res/ind].
+
+optional arguments:
+
+-h, --help show this help message and exit
+-m similarity_matrix.txt, --matrix similarity_matrix.txt
+NumPy matrix containing existing packing similarity
+results.
+-ns 25, --n_structures 25
+Number of structures to take from input set.
+-nm 15, --n_molecules 15
+Size of molecular packing shell to use for analysis
+(must be consistent with input matrix, if used).
+-o Flag for whether to save packing similarity results
+(text file and mol2 overlays).
+--allow_molecular_differences
+Flag for whether to allow for molecular differences
+between structures (e.g. for salts).
+--clustering_type {complete,single,average}
+Type of clustering to employ
+-s, --strip Strip all terminal atoms and alkyl chains, up to any
+hetero atom (O, N, S) or cyclic atom. This cuts the
+
+
+molecule down to core structural features and may
+reveal more general structural similarities, including
+those between molecules with different conformations.
+
+-ct 0.5, --conf_tol 0.
+RMSD threshold for considering two conformations to be
+the same (when merging at level 1).
+-at 25, --angle_tol 25
+Tolerance for angles (in degrees) used by packing
+similarity.
+-dt 0.25, --dist_tol 0.
+Fractional tolerance for distances (0.0 - 1.0) used by
+packing similarity.
+```
+
+## Basic usage (in a command prompt):
+
+```cmd
+python Packing_Similarity_Dendrogram.py input_file
+```
+
+where `input_file` should be a format recognised by the CSD Python API (e.g. cif, res, mol2, gcd or ind). The default
+output of the tool consists of figures (as .png files) of the packing-similarity dendrogram and a heat map of the
+packing similarity between the structures. The file names are prefixed with the stem of the input file ( _i.e._ if
+roy.gcd is the input, then the similarity dendrogram will be called `roy_packing_similarity_tree.png`).
+
+If a large set of structures are inputted, the top _N_ structures can be selected using the `-nm` options.
+
+Using the `-o` option will result in overlays being saved for each comparison (as .mol2 files). The matrix of
+similarities is also saved as a raw numpy matrix ( _i.e._ as
+`roy_similarity_matrix.txt`), and can be read back in (with `-m`) to skip the packing-similarity analysis, if, for
+example, a different clustering algorithm is desired or part of the script has been changed.
+
+The `--allow_molecular_differences` option can be used when comparing crystal structures of closely related moleculese.g. salts and free forms.
+
+The `-s` option will strip all terminal atoms and carbon atom chains up to hetero atoms (e.g.
+the methyl of a methoxy will be removed), which may be useful for identifying more coarse- grained similarity that
+ignores small changes in the periphery of the molecule. A new cif is created containing the stripped molecule and
+analysis is performed using this file. Note that disordered experimental structures present in the input database will
+cause problems when this option is selected and should be removed from the file.
+
+The remaining options control the packing-similarity settings, such as number of molecules, thresholds etc. Matches of
+only one molecule between structures may not necessarily
+
+correspond to good agreement between two conformations (as the distance and angle thresholds can be more coarse-grained
+than an RMSD tolerance). The `-ct` option can therefore be used to specific an RMSD tolerance for whether two
+conformations are matched and in, turn, whether two clusters that have only one molecule in common should be merged at
+that level.
+
+## Understanding the Dendrogram:
+
+An example dendrogram (produced using CBZdataset.ind and related files; see supporting material
+of [Cryst. Growth Des., 2009, 9, 1869 _–_ 1888](https://pubs.acs.org/doi/10.1021/cg801056c) ) is shown below in Figure 1.
+To produce this dendrogram, the script compares all 50 structures in the input to each other. For a default cluster
+size of 15, the match can range from 0, which indicates a completely different conformation and no similarity in
+packing, to 15, where the two structures are isostructural.
+
+In comparing the structures different linkage criteria can be used. The _default behaviour_ is to link structures based
+on the smallest difference between them (referred to as single- linkage clustering). Structures are initially grouped
+together based on their best matches, such that if A matches B with 15/15 and C with 14/15 and B matches C with 15/15,
+then A, B and C are all grouped together at level 15. Groups are then merged with the group which they share the best
+overlay, based on a single structure in each group. This means that if D has an 8/15 molecule match with B, then the
+group A, B, C will join with D at level 8, even if A and C had smaller similarities with D.
+
+At the highest level ( _e.g._ 15 in the example below) isostructural crystal structures are listed together ( _e.g._
+“CBZDBF18, CBZDBF19”). At lower levels, where groups or structures share a similarity then the parent groups are draw
+one level higher (with blue dots) and then join at the level they match. For example, structures 37 and 4 at the bottom
+of the plot have 5 / molecules in common, so both are drawn at level 6 and then merge at level 5. This group of 37 and 4
+then shares a similarity with 32, 09 and 30, as well as all other structures apart from 31 at level 3, so all of them
+merge.
+
+Alternative clustering modes can be invoked by using `--clustering_type`, which takes one of three arguments: single,
+complete and average. Single linkage clustering is the default discussed above. Complete clustering will instead only
+join two clusters at the level where all the structures in the two clusters match, which is equivalent to joining
+clusters at the worst or lowest level connecting them. Figure 2 shows the same database as Figure 1 using this approach,
+which gives more definite clustering but will also hide some similarities due to the hierarchical nature of the
+clustering.
+
+
+**Figure 1** : Example dendrogram based on a database of carbamazepine solid forms using the default single-linkage
+clustering.
+
+
+**Figure 2** : Example dendrogram based on a database of carbamazepine solid forms using complete- linkage clustering.
+
+The final option for clustering type is to merge clusters based on their average packing- similarity agreement. The
+result for carbamazepine is shown in Figure 3. With average linkage clustering, two cluster can join at non-integer
+values and therefore the plot may become more complex visually.
+
+There are several alternative schemes for linking clusters ( _e.g._ weighting the average based on the cluster sizes)
+that could also be implemented.
+
+
+**Figure 3** : Example dendrogram based on a database of carbamazepine solid forms linking clusters based on their
+average packing similarity.
+
+## Caveats:
+
+The tool considers only the heaviest component in each crystal structure to enable comparison of multi-component forms
+with pure forms. Results for _Z_ ’ > 1 systems may not fully reflect differences between structures as the best match is
+retained by default in the packing-similarity analysis.
+
+## Authors:
+_Anthony Reilly(2016)_
+
+> For feedback or to report any issues please contact [support@ccdc.cam.ac.uk](support@ccdc.cam.ac.uk)
\ No newline at end of file
diff --git a/scripts/packing_similarity_dendrogram/assets/dendogram_figure_1.png b/scripts/packing_similarity_dendrogram/assets/dendogram_figure_1.png
new file mode 100644
index 0000000..e665020
Binary files /dev/null and b/scripts/packing_similarity_dendrogram/assets/dendogram_figure_1.png differ
diff --git a/scripts/packing_similarity_dendrogram/assets/dendogram_figure_2.png b/scripts/packing_similarity_dendrogram/assets/dendogram_figure_2.png
new file mode 100644
index 0000000..5e9d899
Binary files /dev/null and b/scripts/packing_similarity_dendrogram/assets/dendogram_figure_2.png differ
diff --git a/scripts/packing_similarity_dendrogram/assets/dendogram_figure_3.png b/scripts/packing_similarity_dendrogram/assets/dendogram_figure_3.png
new file mode 100644
index 0000000..d3fbc96
Binary files /dev/null and b/scripts/packing_similarity_dendrogram/assets/dendogram_figure_3.png differ
diff --git a/scripts/packing_similarity_dendrogram/packing_similarity_dendogram.py b/scripts/packing_similarity_dendrogram/packing_similarity_dendogram.py
new file mode 100644
index 0000000..4356fd3
--- /dev/null
+++ b/scripts/packing_similarity_dendrogram/packing_similarity_dendogram.py
@@ -0,0 +1,457 @@
+#
+# 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-22: created by Anthony Reilly, The Cambridge Crystallographic Data Centre
+# 2016-12-06: updated by Anthony Reilly, The Cambridge Crystallographic Data Centre
+#
+
+"""Packing_Similarity_Dendrogram.py - Construct a dendrogram for an input set of structures based on packing-similarity
+analysis
+"""
+
+import sys
+import argparse
+import matplotlib
+
+matplotlib.use('Agg')
+from ccdc.io import EntryReader, CrystalWriter, MoleculeWriter
+from ccdc.crystal import PackingSimilarity
+import numpy as np
+
+import matplotlib.pyplot as plt
+import os
+
+
+def strip_terminal(name, reader):
+ """
+ Iteratively removes terminal atoms and chains from molecules, stopping at atoms in the allowed list or in rings
+ Saves the resulting structures to "stripped.cif" for use in code
+ """
+ allowed_terminal_atoms = ["S", "O", "N"]
+ xtal_writer = CrystalWriter(name + "_stripped.cif")
+
+ for entry in reader:
+ crystal = entry.crystal
+ crystal.assign_bonds()
+ molecule = crystal.molecule
+ # Zealously strip atoms with only one bond, unless allowed by filter above (e.g. C=O)
+ atom_diff = -1
+ while atom_diff != 0:
+ n_atoms_before = len(molecule.atoms)
+
+ for atom in molecule.atoms:
+ if len(atom.bonds) == 1 and atom.atomic_symbol not in allowed_terminal_atoms:
+ molecule.remove_atom(atom)
+ n_atoms_after = len(molecule.atoms)
+ atom_diff = n_atoms_after - n_atoms_before
+
+ crystal.molecule = molecule
+ xtal_writer.write(crystal)
+
+
+# Cluster functions
+def compare_clusters(c1, c2, relations, mode):
+ # Set sensible bounds on the cluster levels
+ if mode == 'complete':
+ level = 10000
+ for id1 in c1['identifiers']:
+ for id2 in c2['identifiers']:
+ l = relations[id1][id2]
+ level = min(level, l)
+ return level
+ elif mode == 'single':
+ level = -1
+ for id1 in c1['identifiers']:
+ for id2 in c2['identifiers']:
+ l = relations[id1][id2]
+ level = max(level, l)
+ return level
+ elif mode == 'average':
+ l = []
+ for id1 in c1['identifiers']:
+ for id2 in c2['identifiers']:
+ l.append(relations[id1][id2])
+ level = sum(l) / len(l)
+ return level
+
+
+def merge_equal_levels(cluster):
+ children = cluster['children']
+ for child in children:
+ merge_equal_levels(child)
+ pos = 0
+ while pos < len(children):
+ child = children[pos]
+ if child['level'] == cluster['level']:
+ children.remove(child)
+ children.extend(child['children'])
+ else:
+ pos += 1
+
+
+def merge_clusters(c1, c2, level):
+ new_cluster = {'level': level, 'identifiers': c1['identifiers'] + c2['identifiers'], 'children': [c1, c2]}
+ return new_cluster
+
+
+def plot_dendrogram(cluster_list, n_ps_mols, filename):
+ """
+ Function for producing a dendrogram from an input cluster hierarchy
+ """
+
+ def get_terminal_count(cluster, total):
+ """
+ Find the number of terminals in the current cluster
+ """
+ if not cluster['children']:
+ total += 1
+ else:
+ for child in cluster['children']:
+ total = get_terminal_count(child, total)
+ return total
+
+ def get_terminals(cluster, terminals, hs):
+ """
+ Find the heights of the cluster terminals
+ """
+ if not cluster['children']:
+ terminals.append(float(hs[str(cluster['identifiers'])]))
+ else:
+ for child in cluster['children']:
+ terminals = get_terminals(child, terminals, hs)
+ return terminals
+
+ def assign_y_positions(cluster, c, terminal_y_positions):
+ """
+ Assign heights to each of the terminal positions in the current cluster
+ """
+ if not cluster['children']:
+ c += 1
+ terminal_y_positions.update({str(cluster['identifiers']): c})
+ else:
+ for child in cluster['children']:
+ terminal_y_positions, c = assign_y_positions(child, c, terminal_y_positions)
+ return terminal_y_positions, c
+
+ def get_midpoint(cluster, hs):
+ """
+ Find the midpoint of the current cluster
+ """
+ terminal_list = []
+ terminal_list = get_terminals(cluster, terminal_list, hs)
+ midpoint = sum(terminal_list) / len(terminal_list)
+ return midpoint
+
+ def plot_tree(cluster, x_start_positions, y_start_positions, hs):
+ """
+ Recursive function to plot a dendrogram
+ """
+ node = cluster['level']
+
+ if cluster['children']:
+ terminal_list = []
+ terminal_list = get_terminals(cluster, terminal_list, hs)
+ midpoint = sum(terminal_list) / len(terminal_list)
+
+ xpos = [x_start_positions[1], node]
+ ypos = [y_start_positions[1], midpoint]
+ else:
+ if len(cluster['identifiers']) != 1:
+ # Starting groups that are merged start at the highest level to indicate this
+ xpos = [x_start_positions[1], node]
+ else:
+ # Point joins at the next level
+ xpos = [x_start_positions[1], x_start_positions[1] + 1]
+ ypos = [y_start_positions[1], hs[str(cluster['identifiers'])]]
+
+ if node != 0:
+ # Plotting x across and y up or down
+ plt.plot(xpos, [ypos[1], ypos[1]], marker=None, linestyle='-', linewidth=1.0,
+ color="Black", zorder=2)
+ plt.plot([xpos[0], xpos[0]], ypos, marker=None, linestyle='-', linewidth=1.0,
+ color="Black", zorder=2)
+
+ if cluster['children']:
+ # Now repeat for each child
+ for child in cluster['children']:
+ plot_tree(child, xpos, ypos, hs)
+ else:
+ # Plot terminal point and add structure indices
+ plt.scatter(xpos[1], ypos[1], color="Blue", zorder=3, s=10)
+ plt.annotate(str(",".join(cluster['identifiers'])), xy=(xpos[1], ypos[1]),
+ xytext=(xpos[1] + 0.15, ypos[1]), verticalalignment='center', fontsize='3')
+
+ # Setup tree plotting by getting each terminal's height
+ heights, count = assign_y_positions(cluster_list[0], 0, {})
+
+ # Set start of the tree - middle of the plot and 1,1
+ xpositions = [1, 1]
+ ypositions = [1, get_midpoint(cluster_list[0], heights)]
+
+ # Plot tree
+ plot_tree(cluster_list[0], xpositions, ypositions, heights)
+
+ # Plot formatting
+ ax = plt.axes()
+ ax.set_frame_on(False)
+ ax.axes.get_yaxis().set_visible(False)
+ ax.spines['top'].set_visible(False)
+ ax.spines['bottom'].set_visible(False)
+ ax.spines['left'].set_visible(False)
+ ax.spines['right'].set_visible(False)
+
+ levels = range(n_ps_mols, -1, -1)
+ highlighted_levels = range(n_ps_mols, -1, -1)
+ # highlighted_levels.append(0)
+
+ for level in highlighted_levels:
+ plt.plot([level, level], [0, count], "--", linewidth=0.5, color="Gray", zorder=1)
+
+ # Pad the plot to have enough space for structure indices
+ plt.xlim(-1, n_ps_mols + 5)
+ plt.ylim(0, count + 1)
+ ax.set_xticks(levels)
+ ax.tick_params(axis='x', bottom='off', top='off')
+ plt.xlabel('Packing Similarity / ' + str(n_ps_mols) + ' Molecules', fontsize='large')
+ # Save output
+ plt.savefig(filename + "_packing_similarity_tree.png", dpi=1000, bbox_inches='tight')
+ print("Packing tree diagram saved to " + filename + "_packing_similarity_tree.png")
+
+
+def main(input_file, matrix_file, n_ps_mols, output_ps_results, conf_threshold, ps_angles, ps_distances, strip,
+ n_struct, allow_mol_diff, cluster_mode):
+ # Initialise Packing Similarity
+ ps = PackingSimilarity()
+ ps.settings.ignore_hydrogen_positions = True
+ ps.settings.ignore_bond_types = True
+ ps.settings.match_entire_packing_shell = False
+ # Deal with e.g. salt forms
+ if allow_mol_diff:
+ ps.settings.allow_molecular_differences = True
+ ps.settings.ignore_hydrogen_counts = True
+ ps.settings.ignore_bond_counts = True
+ else:
+ ps.settings.allow_molecular_differences = False
+ ps.settings.ignore_hydrogen_counts = False
+ ps.settings.ignore_bond_counts = False
+
+ # Deal with solvates
+ ps.settings.ignore_smallest_components = True
+ ps.settings.packing_shell_size = n_ps_mols
+ ps.settings.angle_tolerance = ps_angles
+ ps.settings.distance_tolerance = ps_distances
+
+ refcodes = []
+
+ input_name = input_file.rsplit(".")[0]
+ print("--------------------------------------------------------")
+
+ if not matrix_file:
+ # Read in the input file
+ print("Reading input database/gcd/structures:", input_file)
+
+ if not strip:
+ structure_reader = EntryReader(input_file)
+ else:
+ structure_reader = EntryReader(input_file)
+ strip_terminal(input_name, structure_reader)
+ structure_reader.close()
+ structure_reader = EntryReader(input_name + "_stripped.cif")
+
+ if n_struct:
+ # noinspection PyTypeChecker
+ structure_size = min(len(structure_reader), n_struct)
+ else:
+ # noinspection PyTypeChecker
+ structure_size = len(structure_reader)
+
+ # Initialise matrix
+ matrix = np.zeros((structure_size, structure_size))
+
+ print("Generating matrix of packing similarities")
+
+ overlay_folder = input_name + "_overlays"
+ if output_ps_results:
+ g = open("packing_similarity_results.txt", "w")
+ g.write("Packing Similarity Analysis for: " + input_file + "\n")
+ if not os.path.exists(overlay_folder):
+ os.makedirs(overlay_folder)
+
+ for i in range(0, structure_size):
+ refcodes.append(str(structure_reader[i].identifier))
+
+ for i in range(0, structure_size):
+ entry_i = structure_reader[i]
+ crystal_i = entry_i.crystal
+ refcodes.append(str(structure_reader[i].identifier))
+
+ for j in range(i, structure_size):
+ if i == j:
+ matrix[i, i] = ps.settings.packing_shell_size
+ continue
+ entry_j = structure_reader[j]
+ crystal_j = entry_j.crystal
+
+ result = ps.compare(crystal_i, crystal_j)
+
+ if result is not None:
+ if result.nmatched_molecules != 1:
+ matrix[i, j] = result.nmatched_molecules
+ matrix[j, i] = result.nmatched_molecules
+ # For single-molecule matches enforce a threshold for conformation RMSDs
+ elif result.rmsd < conf_threshold:
+ matrix[i, j] = result.nmatched_molecules
+ matrix[j, i] = result.nmatched_molecules
+ else:
+ matrix[i, j] = 0
+ matrix[j, i] = 0
+ if output_ps_results:
+ g.write(entry_i.identifier + " " + entry_j.identifier + ": " + str(int(matrix[i, j])) +
+ " molecules \n")
+
+ overlay_writer = MoleculeWriter(os.path.join(overlay_folder,
+ "overlay_{0}_{1}.mol2".format(refcodes[i],
+ refcodes[j])))
+ mols = result.overlay_molecules()
+ for mol in mols:
+ overlay_writer.write(mol)
+
+ else:
+ matrix[i, j] = 0
+ matrix[j, i] = 0
+ if output_ps_results:
+ g.write(entry_i.identifier + " " + entry_j.identifier + ": " + str(int(matrix[i, j])) +
+ " molecules \n")
+
+ np.savetxt(input_name + '_similarity_matrix.txt', matrix, delimiter=',')
+ print("Packing similarity matrix saved to similarity_matrix.txt")
+ if output_ps_results:
+ g.close()
+ else:
+ # Use the input matrix
+ print("Reading input matrix:", matrix_file)
+ matrix = np.loadtxt(matrix_file, delimiter=',')
+ # Check the matrix is consistent with options specified and the input structures
+ if int(matrix[1, 1]) != n_ps_mols:
+ print("Error - input matrix and requested packing-shell size do not match")
+ sys.exit(1)
+ structure_reader = EntryReader(input_file)
+ # noinspection PyTypeChecker
+ structure_size = len(structure_reader)
+ for i in range(0, structure_size):
+ refcodes.append(str(structure_reader[i].identifier))
+
+ structure_reader.close()
+ if len(matrix) != structure_size:
+ print("Error - input matrix does not contain the same number of structures as inputted")
+ sys.exit(1)
+
+ print("--------------------------------------------------------")
+
+ # Set structure label based on location in input file
+ labels = []
+ for i in range(0, structure_size):
+ labels.append(str(i + 1))
+
+ # Generate a look-up table of results from matrix
+ relations = {}
+ for i in range(0, structure_size):
+ relations.update(
+ {refcodes[i]: {refcodes[x]: matrix[i, x].tolist() for x in range(0, structure_size)}})
+
+ # Generate a cluster hierarchy - populate it initially with every structure
+ cluster_list = []
+ for i in range(0, structure_size):
+ cluster_list.append({'level': n_ps_mols, 'identifiers': [refcodes[i]], 'children': []})
+
+ # Merge the structures - getting best match for each structure
+ while len(cluster_list) > 1:
+ best_score = -1
+ best_cluster1 = None
+ best_cluster2 = None
+ for i in cluster_list:
+ for j in cluster_list:
+ if i == j:
+ continue
+ score = compare_clusters(i, j, relations, cluster_mode)
+ if score > best_score:
+ best_score = score
+ best_cluster1 = i
+ best_cluster2 = j
+ new_cluster = merge_clusters(best_cluster1, best_cluster2, best_score)
+ cluster_list.remove(best_cluster1)
+ cluster_list.remove(best_cluster2)
+ cluster_list.append(new_cluster)
+
+ # Tidy cluster hierarchy by merging groups with equal matches
+ merge_equal_levels(cluster_list[0])
+
+ # Plot a heat map of the PS matrix
+ x = np.arange(0, structure_size + 1, 1)
+ y = np.arange(0, structure_size + 1, 1)
+
+ plot = plt.pcolor(x, y, matrix, cmap=plt.get_cmap('rainbow', (n_ps_mols - 1)), vmin=1, vmax=n_ps_mols)
+ plt.xticks(np.arange(0, structure_size + 1, 5) - 0.5, np.arange(0, structure_size + 1, 5))
+ plt.yticks(np.arange(0, structure_size + 1, 5) - 0.5, np.arange(0, structure_size + 1, 5))
+ cb = plt.colorbar(plot, ticks=range(1, 16))
+
+ plt.xlim(0, structure_size)
+ plt.ylim(0, structure_size)
+ ax = plt.gca()
+ ax.set_xlabel("Structure Index", fontsize='x-large')
+ ax.set_ylabel("Structure Index", fontsize='x-large')
+ cb.set_label('Packing Similarity /' + str(n_ps_mols) + ' Molecules', fontsize='x-large')
+ plt.savefig(input_name + "_heat_map.png", dpi=300)
+ print("Packing similarity heat map saved to " + input_name + "_heat_map.png")
+ plt.close()
+
+ # Plot a dendrogram
+ plot_dendrogram(cluster_list, n_ps_mols, input_name)
+
+ print("--------------------------------------------------------")
+
+ sys.exit()
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter, description=__doc__)
+ parser.add_argument('input_file', help='Set of structures to perform analysis on [.mol2/cif/res/ind]')
+ parser.add_argument('-m', '--matrix', type=str, help='NumPy matrix containing existing packing similarity results.',
+ metavar='similarity_matrix.txt')
+ parser.add_argument('-ns', '--n_structures', type=int, help='Number of structures to take from input set.',
+ metavar='25')
+ parser.add_argument('-nm', '--n_molecules', type=int, default=15,
+ help='Size of molecular packing shell to use for analysis '
+ '(must be consistent with input matrix, if used).', metavar="20")
+ parser.add_argument('-o', action="store_true",
+ help='Flag for whether to save packing similarity results (text file and mol2 overlays).')
+ parser.add_argument('--allow_molecular_differences', action="store_true",
+ help='Flag for whether to allow for molecular differences between structures (e.g. for salts).')
+ parser.add_argument('--clustering_type', choices=['complete', 'single', 'average'], default='single',
+ help='Type of clustering to employ')
+ parser.add_argument('-s', '--strip', action="store_true",
+ help='Strip all terminal atoms and alkyl chains, up to any hetero atom (O, N, S) or cyclic '
+ 'atom. This cuts the molecule down to core structural features and may reveal'
+ ' more general structural similarities, including those between molecules with different '
+ 'conformations.')
+ parser.add_argument('-ct', '--conf_tol', type=float, default=0.5,
+ help='RMSD threshold for considering two conformations to be the same (when merging at level 1)'
+ '.', metavar="0.5")
+ parser.add_argument('-at', '--angle_tol', type=float, default=25, metavar="25",
+ help="Tolerance for angles (in degrees) used by packing similarity.")
+ parser.add_argument('-dt', '--dist_tol', type=float, default=0.25, metavar="0.25",
+ help="Fractional tolerance for distances (0.0 - 1.0) used by packing similarity.")
+ args = parser.parse_args()
+ if not os.path.isfile(args.input_file):
+ parser.error('%s not found.' % args.input_file)
+ if args.matrix:
+ if not os.path.isfile(args.matrix):
+ parser.error('%s not found.' % args.matrix)
+
+ main(args.input_file, args.matrix, args.n_molecules, args.o, args.conf_tol, args.angle_tol,
+ args.dist_tol, args.strip, args.n_structures, args.allow_molecular_differences,
+ args.clustering_type)