Merge pull request #69 from ccdc-opensource/surface_charge_script

Added Topology Surface Charge SCI-1404

Author: Alex-AMC

scripts/surface_charge/ReadMe.md

@@ -11,7 +11,7 @@ Charges are currently calculated using the Gasteiger charge model. Further devel
 
 Example Output:
 
-![Example Output](assets/example_output.png)
+![Example Output](assets/example_output_hxacan28.png)
 
 > **Note** - When comparing charges for non-CSD structures and structures from mol2 files the values might be different as the bonding might not be the same. When importing a mol2 file the bonding and charges may have to be calculated on the fly, whereas this information is assigned for CSD entries.
 

scripts/surface_charge/assets/example_output.png

@@ -8,10 +8,9 @@
 # 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
+from typing import List, Tuple, Dict
 
 import numpy as np
 from ccdc.io import CrystalReader
@@ -29,30 +28,99 @@ def __init__(self, crystal, use_existing_charges: bool = False):
         self.crystal = crystal
         self.surface = None
         self._surface_charge = None
+        self.surface_atom_charge = np.nan
+        self.node_projected_charge = np.nan
+        self.node_representative_charge = np.nan
+        self.triangles_properties = dict()
+
+    @staticmethod
+    def sum_atom_charge(atoms: List[object]) -> float:
+        return np.round(np.sum([atom.partial_charge for atom in atoms]), 3)
+
+    def get_node_charge(self):
+        self.node_charge_dictionary = {}
+        node_list = list(self.surface.topology.nodes)
+        for node, atoms in self.surface.surface_node_atom_contacts.items():
+            node_index = node_list.index(node)
+            total_node_charge = 0
+            if len(atoms) > 0:
+                total_node_charge = self.sum_atom_charge(atoms)
+            self.node_charge_dictionary[node_index] = total_node_charge
+
+    def calculate_triangles_properties(self,
+                                       tri_index: List[Tuple[int, int, int]]) -> None:
+        surface_area = self.surface.descriptors.surface_area
+        self.triangles_properties = {}
+        triangle_areas = self.calculate_area_of_triangles(list(self.surface.topology.triangles))
+
+        for node_index, triangle_area in zip(tri_index, triangle_areas):
+            average_triangle_charge = np.mean([self.node_charge_dictionary[i] for i in node_index])
+            triangle_representation = triangle_area / surface_area
+            projected_charge = 0
+            if np.isclose(triangle_area, 0.0) is False:
+                projected_charge = average_triangle_charge / triangle_area
+            self.triangles_properties[tuple(node_index)] = {'Average Charge': average_triangle_charge,
+                                                            'Triangle Area': triangle_area,
+                                                            'Percentage Area': triangle_representation,
+                                                            'Node Representative Charge': average_triangle_charge * triangle_representation,
+                                                            'Node Projected Surface Charge': projected_charge}
+
+    def calculate_node_charges(self):
+        tri_index = self.calculated_node_index_values(list(self.surface.topology.nodes),
+                                                      list(self.surface.topology.triangles))
+        self.get_node_charge()
+        self.calculate_triangles_properties(tri_index)
+        self.representative_charge = np.sum(
+            [triangle['Node Representative Charge'] for triangle in self.triangles_properties.values()])
+        self.node_charges = np.sum([triangle['Average Charge'] for triangle in self.triangles_properties.values()])
+        return self.representative_charge
+
+    @staticmethod
+    def calculate_length(origin: np.ndarray, target: np.ndarray) -> float:
+        """Returns distance between target and origin"""
+        if not isinstance(origin, np.ndarray) or not isinstance(target, np.ndarray):
+            raise TypeError("Please supply numpy arrays for lengths.")
+        return np.linalg.norm(target - origin)
+
+    @staticmethod
+    def compute_simplex_area(simplex: np.ndarray) -> float:
+        vec_1 = simplex[1] - simplex[0]
+        vec_2 = simplex[2] - simplex[0]
+        return np.linalg.norm(np.cross(vec_1, vec_2)) / 2
+
+    def calculate_area_of_triangles(self, triangles: List) -> List:
+        """ Calculates area of individual triangles from node positions using Heron's formula"""
+        return [self.compute_simplex_area(np.array(triangle)) for triangle in triangles]
+
+    @staticmethod
+    def calculated_node_index_values(nodes: List, triangles: List) -> List:
+        """
+         Calculate index of all triangles for nodes
+
+        :param nodes: list of nodes [x,y,z]
+        :param triangles: list of triangles that contain nodes index values
+        """
+        search_dictionary = {e: i for i, e in enumerate(nodes)}
+        return [[search_dictionary[n] for n in tri] for tri in triangles]
 
     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)
+            self.surface_atom_charge = self.sum_atom_charge(atoms=self.surface.surface_atoms)
+            self.node_representative_charge = self.calculate_node_charges()
             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_node_charges = []
+        self.surface_areas = []
+        self.surface_node_representative_charge = []
+        self.surface_atom_charges = []
         self.surface_charges_per_area = []
         self.surface_charges = []
         self.projected_area = []
@@ -84,9 +152,12 @@ 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.surface_atom_charges.append(charges.surface_atom_charge)
+            self.surface_node_charges.append(charges.node_charges)
+            self.surface_node_representative_charge.append(charges.node_representative_charge)
+
             self.projected_area.append(charges.surface.descriptors.projected_area)
-            self.surface_charges_per_area.append(charges.surface_charge_per_area)
+            self.surface_areas.append(charges.surface.descriptors.surface_area)
 
     def make_report(self):
         html_content = self.generate_html_table()
@@ -114,11 +185,13 @@ def generate_html_table(self):
             <h2>Calculation Results</h2>
             <table>
                 <tr>
-                    <th>hkl</th>
+                    <th width="80px">hkl</th>
                     <th>Offset</th>
-                    <th>Projected Area</th>
-                    <th>Surface Charge*</th>
-                    <th>Surface Charge per Projected Area</th>
+                    <th>P<sub>Area</sub>-Projected Area &#8491;<sup>2</sup></th>
+                    <th>S<sub>Area</sub>-Surface Area &#8491;<sup>2</sup></th>
+                    <th>Total Atom Surface Charge</th>
+                    <th>Total Atom Surface Charge/P<sub>A</sub></th>
+                    <th>Topological Surface Charge/ S<sub>Area</sub></th>
                 </tr>
         """
 
@@ -130,15 +203,18 @@ def generate_html_table(self):
                     <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>
+                    <td>{self.surface_areas[i]:.3f}</td>
+                    <td>{self.surface_atom_charges[i]:.3f}</td>
+                    <td>{self.surface_atom_charges[i] / self.projected_area[i]:.4f}</td>
+                    <td>{self.surface_node_representative_charge[i]:.4f}</td>
                 </tr>
             """
 
         # HTML Table Footer
         html += """
             </table>
             <p><i> *-Surface charge is based on gasteiger partial charges <a href="https://www.sciencedirect.com/science/article/pii/S0040403901949779?via%3Dihub">10.1016/S0040-4039(01)94977-9</a></i> </p>
+            <p> Topological surface charge is defined as the average triangle charge on the surface multiplied by the % area contribution towards the total surface. </p> 
         </body>
         </html>
         """