diff --git a/scripts/surface_charge/ReadMe.md b/scripts/surface_charge/ReadMe.md index 5f4c140..b186a20 100644 --- a/scripts/surface_charge/ReadMe.md +++ b/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. diff --git a/scripts/surface_charge/assets/example_output.png b/scripts/surface_charge/assets/example_output.png deleted file mode 100644 index 3f9205a..0000000 Binary files a/scripts/surface_charge/assets/example_output.png and /dev/null differ diff --git a/scripts/surface_charge/assets/example_output_hxacan28.png b/scripts/surface_charge/assets/example_output_hxacan28.png new file mode 100644 index 0000000..09c0505 Binary files /dev/null and b/scripts/surface_charge/assets/example_output_hxacan28.png differ diff --git a/scripts/surface_charge/surface_charge_calculator.py b/scripts/surface_charge/surface_charge_calculator.py index 498c891..80cbbc7 100644 --- a/scripts/surface_charge/surface_charge_calculator.py +++ b/scripts/surface_charge/surface_charge_calculator.py @@ -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):

Calculation Results

- + - - - + + + + + """ @@ -130,8 +203,10 @@ def generate_html_table(self): - - + + + + """ @@ -139,6 +214,7 @@ def generate_html_table(self): html += """
hklhkl OffsetProjected AreaSurface Charge*Surface Charge per Projected AreaPArea-Projected Area Å2SArea-Surface Area Å2Total Atom Surface ChargeTotal Atom Surface Charge/PATopological Surface Charge/ SArea
{hkl} {offset:.2f} {self.projected_area[i]:.3f}{self.surface_charges[i]:.3f}{self.surface_charges_per_area[i]:.4f}{self.surface_areas[i]:.3f}{self.surface_atom_charges[i]:.3f}{self.surface_atom_charges[i] / self.projected_area[i]:.4f}{self.surface_node_representative_charge[i]:.4f}

*-Surface charge is based on gasteiger partial charges 10.1016/S0040-4039(01)94977-9

+

Topological surface charge is defined as the average triangle charge on the surface multiplied by the % area contribution towards the total surface.

"""