# # 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-04: created by Peter A. Wood, the Cambridge Crystallographic Data Centre # 2015-07-24: made available by the Cambridge Crystallographic Data Centre # import os import matplotlib.pyplot as plot from ccdc import search from ccdc.diagram import DiagramGenerator from mercury_interface import MercuryInterface from ccdc import conformer from PIL import Image, ImageDraw, ImageOps def make_diagram(mol): # 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 try: img = molecule_diagram_generator.image(mol) except RuntimeError: img = make_blank_image('No 2D Diagram\nAvailable') fname = os.path.join(output_dir, '%s_diagram.png' % crystal.identifier) if img: img.save(fname) return os.path.basename(fname) def make_blank_image(text=None, size=200, border=2): # Generate a blank image (for 2D diagram failure) img = Image.new('RGB', (size-border, size-border), 'white') draw = ImageDraw.Draw(img) draw.text((10,10), text, fill='black') img = ImageOps.expand(img, border=border, fill='black') return img def write_histogram(hist_name, data, ref_val, settings): # make file name file_name = str(os.path.join(output_dir, '%s_histogram.png' % hist_name)) # Work out where the histogram bins should be placed on the histogram xs = [settings["low"] + i * settings["bin_width"] for i in range(len(data))] # Set up the basic plot fig = plot.figure() ax = fig.add_subplot(111) ax.bar(xs, data, settings["bin_width"], color='g') ax.set_xlim(settings["low"], settings["high"]) ax.set_xlabel(settings["xlabel"]) ax.set_ylabel(settings["ylabel"]) ax.set_title(settings["title"]) # Add the reference value line ax.add_line(plot.Line2D((ref_val, ref_val), plot.ylim(), color='r')) # Add the text describing the query value next to the reference line. ref_val_str = 'Value in query: %.3f' % ref_val middle = (settings["low"] + settings["high"]) / 2.0 ymin, ymax = plot.ylim() xoffset = (settings["high"] - settings["low"]) / 100 yoffset = (ymin - ymax) / 80 ref_y = ymax + yoffset if ref_val > middle: ref_x = ref_val - xoffset ha = 'right' else: ref_x = ref_val + xoffset ha = 'left' ax.text(ref_x, ref_y, ref_val_str, color='r', ha=ha, va='top') plot.savefig(file_name) plot.close(fig) def write_bond_distance_histogram(bond, bond_name): # Write bond distance histogram to a file settings = {} settings["low"] = min(bond.value, bond.minimum) - 0.1 settings["high"] = max(bond.value, bond.maximum) + 0.1 diff = settings["high"] - settings["low"] settings["bin_width"] = diff / 40.0 settings["title"] = 'Bond Lengths for %s (Click to Return to Table)' % bond.fragment_label settings["xlabel"] = 'Bond Length' settings["ylabel"] = 'Number of Hits' data = bond.histogram(settings["bin_width"], settings["low"], settings["high"]) write_histogram(bond_name, data, bond.value, settings) def write_valence_angle_histogram(angle, angle_name): # Write valence angle histogram to a file settings = {} settings["low"] = min(angle.value, angle.minimum) - 1 settings["high"] = max(angle.value, angle.maximum) + 1 diff = settings["high"] - settings["low"] settings["bin_width"] = diff / 40.0 settings["title"] = 'Valence Angles for %s (Click to Return to Table)' % angle.fragment_label settings["xlabel"] = 'Valence Angle' settings["ylabel"] = 'Number of Hits' data = angle.histogram(settings["bin_width"], settings["low"], settings["high"]) write_histogram(angle_name, data, angle.value, settings) def write_torsion_angle_histogram(torsion, torsion_name): # Write torsion angle histogram to a file settings = {} settings["low"] = 0 settings["high"] = 180 settings["bin_width"] = 5.0 settings["title"] = 'Torsion Angles for %s (Click to Return to Table)' % torsion.fragment_label settings["xlabel"] = 'Torsion Angle' settings["ylabel"] = 'Number of Hits' data = torsion.histogram(settings["bin_width"], settings["low"], settings["high"]) write_histogram(torsion_name, data, abs(torsion.value), settings) def geometry_analyser(): # Sets up setting for Geometry Analyser engine = conformer.GeometryAnalyser() engine.settings.bond.analyse = True engine.settings.angle.analyse = True engine.settings.torsion.analyse = True engine.settings.ring.analyse = False return engine if __name__ == '__main__': helper = MercuryInterface() entry_id = helper.identifier # Read crystal relating to entry_id from any supported file format entry = helper.current_entry crystal = entry.crystal molecule = crystal.molecule molecule.assign_bond_types(which='unknown') engine = geometry_analyser() tableno = 1 # Define HTML output file and open it html_file = html_file = helper.output_html_file f = open(html_file, "w") # Extract output directory from input html_file string output_dir = os.path.dirname(html_file) good = ' class="good"' med = ' class="med"' bad = ' class="bad"' # Write report header and title f.write('\n'.join([ '' '' '' '' '' '
' '' '' % entry_id, '' % helper.get_ccdc_logo, '
Would You Publish This? (%s)
' '
' ])) # Write diagram and simple crystal part of HTML output file f.write('\n'.join([ '
' '

' % make_diagram(molecule), '
', ])) # Geometry analysis on molecule - showing only unusual bits and bobs checked_mol = engine.analyse_molecule(molecule) # Generating bond length data table unusual_bonds = [] unusual_angles = [] unusual_torsions = [] if len(checked_mol.analysed_bonds) > 0: f.write('

Intramolecular Geometry Check

') f.write('') f.write('') f.write('' '') f.write('') for bond in checked_mol.analysed_bonds: atom1, atom2 = bond.atom_labels bond_name = "Bond_%s-%s" % (atom1, atom2) bond_value = str(bond.value) if bond.unusual: write_bond_distance_histogram(bond, bond_name) file_name = str(os.path.join(output_dir, '%s_histogram.png' % bond_name)) unusual_bonds.append(bond_name) b_z_score = bond.z_score if b_z_score <= 3.0: bond_decision = med else: bond_decision = bad f.write('' % (bond_name, file_name)) f.write( '%.3f' % (atom1, atom2, '-', '-', bond.value, 'z-score', bond_decision, b_z_score)) for angle in checked_mol.analysed_angles: atom1, atom2, atom3 = angle.atom_labels angle_name = "Angle_%s-%s-%s" % (atom1, atom2, atom3) angle_value = str(angle.value) if angle.unusual: write_valence_angle_histogram(angle, angle_name) file_name = str(os.path.join(output_dir, '%s_histogram.png' % angle_name)) unusual_angles.append(angle_name) a_z_score = angle.z_score if a_z_score <= 3.0: angle_decision = med else: angle_decision = bad f.write('' % (angle_name, file_name)) f.write('%.3f' % (atom1, atom2, atom3, '-', angle.value, 'z-score', angle_decision, a_z_score)) for torsion in checked_mol.analysed_torsions: atom1, atom2, atom3, atom4 = torsion.atom_labels torsion_name = "Torsion_%s-%s-%s-%s" % (atom1, atom2, atom3, atom4) torsion_value = str(torsion.value) if torsion.unusual: write_torsion_angle_histogram(torsion, torsion_name) file_name = str(os.path.join(output_dir, '%s_histogram.png' % torsion_name)) unusual_torsions.append(torsion_name) l_den = torsion.local_density if 3.5 <= l_den <= 5.0: torsion_decision = med else: torsion_decision = bad f.write('' % (torsion_name, file_name)) f.write('%.3f' % (atom1, atom2, atom3, atom4, torsion.value, 'local density', torsion_decision, l_den)) if (len(unusual_bonds) + len(unusual_angles) + len(unusual_torsions)) == 0: f.write('') f.write('
TypeAtom 1Atom 2Atom 3Atom 4ValueMetricMetric Value
bond%s%s%s%s%.3f Å%s
angle%s%s%s%s%.3f°%s
torsion%s%s%s%s%.3f°%s
No unusual intramolecular geometric parameters
') crystal.assign_bonds() # Define H-bond donor, halogen-bond donor, acceptor & ring substructures donor_substructure = search.SMARTSSubstructure("[N,O]-H") hal_donor_substructure = search.SMARTSSubstructure("[#6]-[#35,#53,#85]") acceptor_substructure = search.SMARTSSubstructure("[N,n,O,o,S,Cl,#35,I]") phenyl_substructure = search.SMARTSSubstructure("c1ccccc1") # Find and capture all intermolecular H-bonding details inter_hbond_search = search.SubstructureSearch() inter_d_id = inter_hbond_search.add_substructure(donor_substructure) inter_a_id = inter_hbond_search.add_substructure(acceptor_substructure) inter_hbond_search.add_distance_constraint('HA', inter_d_id, 1, inter_a_id, 0, (-5.00, 0.00), 'Intermolecular', vdw_corrected=True) inter_hbond_search.add_distance_measurement('DA', inter_d_id, 0, inter_a_id, 0) inter_hbond_search.add_distance_measurement('DH', inter_d_id, 0, inter_d_id, 1) inter_hbond_search.add_angle_constraint('DHA', inter_d_id, 0, inter_d_id, 1, inter_a_id, 0, (120.0, 180.0)) inter_hbond_hits = inter_hbond_search.search(crystal) # Find and capture all intramolecular H-bonding details intra_hbond_search = search.SubstructureSearch() intra_d_id = intra_hbond_search.add_substructure(donor_substructure) intra_a_id = intra_hbond_search.add_substructure(acceptor_substructure) intra_hbond_search.add_distance_constraint('HA', intra_d_id, 1, intra_a_id, 0, (-5.00, 0.00), 'Intramolecular', vdw_corrected=True) intra_hbond_search.add_distance_measurement('DA', intra_d_id, 0, intra_a_id, 0) intra_hbond_search.add_distance_measurement('DH', intra_d_id, 0, intra_d_id, 1) intra_hbond_search.add_angle_constraint('DHA', intra_d_id, 0, intra_d_id, 1, intra_a_id, 0, (120.0, 180.0)) intra_hbond_hits = intra_hbond_search.search(crystal) # Write H-bonding part of HTML output file f.write('
') f.write('

Intermolecular Interaction Check †

') num_hbonds = len(inter_hbond_hits) + len(intra_hbond_hits) if num_hbonds > 0: f.write('') f.write('\n'.join([ '' '' '' '' '' '' '' '' ])) for hbond in inter_hbond_hits: if 150.0 <= hbond.constraints['DHA'] <= 180.0: inter_decision = good elif 140.0 <= hbond.constraints['DHA']: inter_decision = med else: inter_decision = bad f.write('\n'.join([ '' '' % (hbond.match_atoms()[0].label, hbond.match_atoms()[1].label), '' % (hbond.match_atoms()[2].label, hbond.measurements['DH']), '' % (hbond.constraints['HA'], hbond.measurements['DA']), '%.1f' % (inter_decision, hbond.constraints['DHA']) ])) for hbond in intra_hbond_hits: if 150.0 <= hbond.constraints['DHA'] <= 180.0: intra_decision = good elif 140.0 <= hbond.constraints['DHA']: intra_decision = med else: intra_decision = bad f.write('\n'.join([ '' '' % (hbond.match_atoms()[0].label, hbond.match_atoms()[1].label), '' % (hbond.match_atoms()[2].label, hbond.measurements['DH']), '' % (hbond.constraints['HA'], hbond.measurements['DA']), '%.1f' % (intra_decision, hbond.constraints['DHA']) ])) f.write('
TypeDonorAcceptorD-H (Å)H...A (Å)D...A (Å)D-H...A (°)
intramolecular%s-%s%s%.2f%.3f%.3f
intramolecular%s-%s%s%.2f%.3f%.3f
') else: f.write('

No H-bonds found in this structure.

') f.write('
') # Packing characterisation calculations heavy = 0 for atom in molecule.atoms: if atom.atomic_symbol != 'H' and atom.atomic_symbol != 'D': heavy += 1 theoretical_volume = heavy * 18 * crystal.z_value experimental_volume = crystal.cell_volume if (0.9 * theoretical_volume) <= experimental_volume <= (1.1 * theoretical_volume): volume_decision = good elif (0.85 * theoretical_volume) <= experimental_volume <= (1.15 * theoretical_volume): volume_decision = med else: volume_decision = bad packing_coefficient = crystal.packing_coefficient if molecule.is_organic: moltype = "organic" av_pack_coeff = "0.68(4)" if 0.6 <= packing_coefficient <= 0.76: packing_decision = good elif 0.56 <= packing_coefficient <= 0.8: packing_decision = med else: packing_decision = bad elif molecule.is_organometallic: moltype = "organometallic" av_pack_coeff = "0.67(5)" if 0.57 <= packing_coefficient <= 0.77: packing_decision = good elif 0.52 <= packing_coefficient <= 0.82: packing_decision = med else: packing_decision = bad else: moltype = "" av_pack_coeff = "" packing_decision = "" # Define probe radius and grid spacing for void calculation probe_radius = 1.2 grid_spacing = 0.7 void_perc = crystal.void_volume(probe_radius, grid_spacing, mode='contact') void_vol = (void_perc * crystal.cell_volume) / 100 if void_perc == 0.0: void_decision = good elif void_perc <= 5.0: void_decision = med else: void_decision = bad # Write packing characterisation part of HTML output file if volume_decision == bad or packing_decision == bad or void_decision == bad: packing_check_decision = bad elif volume_decision == med or packing_decision == med or void_decision == med: packing_check_decision = med else: packing_check_decision = good f.write('\n'.join([ '
' '
Packing Check
' % packing_check_decision, '' ' %d ų ' % ( volume_decision, theoretical_volume), ' %.2f ų ' % (volume_decision, experimental_volume), ' %.3f ' % (packing_decision, packing_coefficient), '' % (moltype, av_pack_coeff), ' %.2f % ' % (void_decision, void_perc), ' %.2f ų ' % (void_decision, void_vol), '
Estimated volume from 18 ų rule
Experimental volume
Packing coefficient
CSD average packing coefficient for %s molecules %s
Calculated void percentage *
Calculated void volume *
', '
', ])) # Write appendix f.write('
') f.write('

Appendix

') f.write('\n'.join([ '

† Hydrogen bond: defined here as a hydrogen-bond donor group (O or N atom', 'bonded to a hydrogen) interacting with an acceptor group (N, O, S, Cl, Br or I).', 'The H···A distance must be less than the sum of the van der Waals radii and the', 'D-H···A angle must be greater than or equal to 120°.

\n\n', ])) f.write('

* Void calculations performed with a probe radius of %.1f Å /' 'and a grid spacing of %.1f Å.

' % (probe_radius, grid_spacing)) for hist in unusual_bonds: file_name = str(os.path.join(output_dir, '%s_histogram.png' % hist)) f.write('\n'.join([ '
\n' % (file_name, hist, file_name), ])) for hist in unusual_angles: file_name = str(os.path.join(output_dir, '%s_histogram.png' % hist)) f.write('\n'.join([ '
\n' % (file_name, hist, file_name), ])) for hist in unusual_torsions: file_name = str(os.path.join(output_dir, '%s_histogram.png' % hist)) f.write('\n'.join([ '
\n' % (file_name, hist, file_name), ])) f.write('
') '' '' # Close HTML file f.close()