Surface Charge Script

scripts/ReadMe.md

@@ -44,6 +44,10 @@ This folder contains scripts submitted by users or CCDC scientists for anyone to
 
 - Calculates the simulated BFDH particle rugosity weighted by facet area.
 
+## Surface Charge
+
+- Calculates the surface charge for a given structure and surface terminations. Runs both from CMD and Mercury.
+
 ## Tips
 
 A section for top tips in using the repository and GitHub.

scripts/surface_charge/ReadMe.md

@@ -0,0 +1,48 @@
+# Surface Charge Calculator
+
+## Summary
+
+This tool returns the surface charge for a given structure and list of supplied hkl and offsets.
+The script provides a GUI that can be used from Mercury or from the command line.
+
+The output is a HTML file with a table for all the all selected surfaces and their associated charge, projected surface areas and normalised surface charge.
+
+Charges are currently calculated using Gasteiger charges. Further development could be made to use user derived charges. Please let me know if that is something you'd like [[email protected]]([email protected]).
+
+Example Output:
+![Example Output](assets/example_output.png)
+
+> **Note** - When comparing charges for structures out of the CSD and from mol2 files the values might be different as the bonding might not be the same. When importing a mol2 the bonding and charges have to be calculated on the fly. Whereas the CSD structures the bonding is pre-assigned.
+
+## Requirements
+
+- Requires a minimum of CSD 2022.2
+
+## Licensing Requirements
+
+- CSD-Particle Licence
+
+## Instructions on Running
+
+- To Run from command line:
+
+```commandline
+# With an activated environment
+> python surface_charge.py
+```
+
+- To run from mercury:
+Add the folder containing the script to your Python API menu. Mercury -> CSD Python API-> Options -> Add Location. Then just select the `surface_charge.py` script from the drop down menu
+![Adding_Locations](assets/adding_location.png)
+![Selecting Scripts](assets/selecting_script.png)
+
+Running either from the command line or Mercury you will get the same interface allowing you to select a refcode from the CSD or input a mol2 file directly.
+
+Example input :
+![Example Input](assets/example_input.png)
+
+## Author
+
+Alex Moldovan (2024)
+
+> For feedback or to report any issues please contact [[email protected]](mailto:[email protected])

scripts/surface_charge/_surface_charge_calculator.py

@@ -0,0 +1,144 @@
+#
+# 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.
+#
+# The following line states a licence feature that is required to show this script in Mercury and Hermes script menus.
+# Created 18/08/2024 by Alex Moldovan
+# ccdc-licence-features not-applicable
+
+import os
+import warnings
+from typing import List, Tuple
+
+import numpy as np
+from ccdc.io import CrystalReader
+from ccdc.particle import Surface
+from ccdc.utilities import HTMLReport
+
+
+class SurfaceCharge:
+    def __init__(self, crystal, use_existing_charges: bool = False):
+        if use_existing_charges:
+            # if crystal.molecule.assign_partial_charges() is False:
+            #     warnings.warn(f"Gasteiger charges could not be assigned to molecule: {crystal.identifier}",
+            #                   RuntimeWarning)
+            raise NotImplementedError("Use existing charges instead. Current implementation only supports Gasteiger.")
+        self.crystal = crystal
+        self.surface = None
+        self._surface_charge = None
+
+    def calculate_surface_charge(self, hkl: Tuple[int, int, int], offset: float):
+        self.surface = Surface(self.crystal, miller_indices=hkl, offset=offset)
+        if self.surface.slab.assign_partial_charges():
+            self._surface_charge = np.round(np.sum([atom.partial_charge for atom in self.surface.surface_atoms]), 3)
+            return
+        warnings.warn(f"Gasteiger charges could not be assigned to molecule: {self.crystal.identifier}",
+                      RuntimeWarning)
+        self._surface_charge = np.nan
+
+    @property
+    def surface_charge(self):
+        if self._surface_charge is None:
+            raise ValueError("Surface charge calculation not yet calculated, run calculate_surface_charge()")
+        return self._surface_charge
+
+    @property
+    def surface_charge_per_area(self):
+        return self.surface_charge / self.surface.descriptors.projected_area
+
+
+class SurfaceChargeController:
+    def __init__(self, structure: str, hkl_and_offsets: List[Tuple[int, int, int, float]],
+                 output_directory: str = None, use_existing_charges: bool = False):
+        self.surface_charges_per_area = []
+        self.surface_charges = []
+        self.projected_area = []
+        self.crystal = None
+        if output_directory is None:
+            output_directory = os.getcwd()
+        self.output_directory = output_directory
+        self.structure = structure
+        self._surfaces = None
+        self.get_structure()
+        self.identifier = self.crystal.identifier
+        self._surfaces = hkl_and_offsets
+        self.use_existing_charges = use_existing_charges
+
+    def get_structure(self):
+        if "." not in self.structure:
+            self.crystal = CrystalReader('CSD').crystal(self.structure)
+        elif ".mol2" in self.structure:
+            self.crystal = CrystalReader(self.structure)[0]
+        else:
+            raise IOError(" \n ERROR : Please supply refcode or  mol2")
+
+    @property
+    def surfaces(self):
+        if self._surfaces:
+            return self._surfaces
+
+    def calculate_surface_charge(self):
+        for surface in self.surfaces:
+            charges = SurfaceCharge(crystal=self.crystal, use_existing_charges=self.use_existing_charges)
+            charges.calculate_surface_charge(hkl=surface[:3], offset=surface[3])
+            self.surface_charges.append(charges.surface_charge)
+            self.projected_area.append(charges.surface.descriptors.projected_area)
+            self.surface_charges_per_area.append(charges.surface_charge_per_area)
+
+    def make_report(self):
+        html_content = self.generate_html_table()
+        output_file = os.path.join(self.output_directory, self.identifier + "_surface_charge.html")
+        with HTMLReport(file_name=output_file,
+                        ccdc_header=True, append=False, embed_images=False,
+                        report_title=f'Surface Charge Calculations - {self.identifier}') as self.report:
+            self.report.write(html_content)
+
+            print(f"Results saved to {output_file}")
+
+    def generate_html_table(self):
+        # HTML Table Header
+        html = """
+        <html>
+        <head>
+            <title>Calculation Results</title>
+            <style>
+                table { width: 100%; border-collapse: collapse; }
+                th, td { border: 1px solid black; padding: 8px; text-align: center; }
+                th { background-color: #f2f2f2; }
+            </style>
+        </head>
+        <body>
+            <h2>Calculation Results</h2>
+            <table>
+                <tr>
+                    <th>hkl</th>
+                    <th>Offset</th>
+                    <th>Projected Area</th>
+                    <th>Surface Charge</th>
+                    <th>Surface Charge per Projected Area</th>
+                </tr>
+        """
+
+        # HTML Table Rows
+        for i, (h, k, l, offset) in enumerate(self.surfaces):
+            hkl = "{" + f"{h}, {k}, {l}" + "}"
+            html += f"""
+                <tr>
+                    <td>{hkl}</td>
+                    <td>{offset:.2f}</td>
+                    <td>{self.projected_area[i]:.3f}</td>
+                    <td>{self.surface_charges[i]:.3f}</td>
+                    <td>{self.surface_charges_per_area[i]:.4f}</td>
+                </tr>
+            """
+
+        # HTML Table Footer
+        html += """
+            </table>
+        </body>
+        </html>
+        """
+        return html