diff --git a/scripts/ReadMe.md b/scripts/ReadMe.md index 6f34c18..023632f 100644 --- a/scripts/ReadMe.md +++ b/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. diff --git a/scripts/surface_charge/ReadMe.md b/scripts/surface_charge/ReadMe.md new file mode 100644 index 0000000..5f4c140 --- /dev/null +++ b/scripts/surface_charge/ReadMe.md @@ -0,0 +1,50 @@ +# Surface Charge Calculator + +## Summary + +This tool returns the total surface charges for a given structure and list of supplied hkl indices and offsets. +The script provides a GUI that can be used from Mercury or from the command line. + +The output is an HTML file with a table for all the selected surfaces and their associated charges, projected surface areas, and normalised surface charges (surface charge per projected area). + +Charges are currently calculated using the Gasteiger charge model. Further development could be made to use user derived charges. Please let us know if that is of interest: [support@ccdc.cam.ac.uk](support@ccdc.cam.ac.uk). + +Example Output: + + + +> **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. + +## Requirements + +- Requires a minimum of CSD 2022.2 + +## Licensing Requirements + +- CSD-Particle Licence + +## Instructions for use + +- 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 select the `surface_charge.py` script from the drop down menu + + + +Running from either the command line or Mercury will show the same interface allowing you to select a refcode from the CSD or input a mol2 file directly. + +Example Input: + + + +## Author + +Alex Moldovan (2024) + +> For feedback or to report any issues please contact [support@ccdc.cam.ac.uk](mailto:support@ccdc.cam.ac.uk) \ No newline at end of file diff --git a/scripts/surface_charge/assets/adding_location.png b/scripts/surface_charge/assets/adding_location.png new file mode 100644 index 0000000..312335a Binary files /dev/null and b/scripts/surface_charge/assets/adding_location.png differ diff --git a/scripts/surface_charge/assets/csd-python-api-logo.png b/scripts/surface_charge/assets/csd-python-api-logo.png new file mode 100644 index 0000000..c4f1242 Binary files /dev/null and b/scripts/surface_charge/assets/csd-python-api-logo.png differ diff --git a/scripts/surface_charge/assets/example_input.png b/scripts/surface_charge/assets/example_input.png new file mode 100644 index 0000000..4256aa7 Binary files /dev/null and b/scripts/surface_charge/assets/example_input.png differ diff --git a/scripts/surface_charge/assets/example_output.png b/scripts/surface_charge/assets/example_output.png new file mode 100644 index 0000000..3f9205a Binary files /dev/null and b/scripts/surface_charge/assets/example_output.png differ diff --git a/scripts/surface_charge/assets/selecting_script.png b/scripts/surface_charge/assets/selecting_script.png new file mode 100644 index 0000000..6dde040 Binary files /dev/null and b/scripts/surface_charge/assets/selecting_script.png differ diff --git a/scripts/surface_charge/surface_charge.py b/scripts/surface_charge/surface_charge.py new file mode 100644 index 0000000..c76ccb1 --- /dev/null +++ b/scripts/surface_charge/surface_charge.py @@ -0,0 +1,270 @@ +# +# 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 (https://orcid.org/0000-0003-2776-3879) + + +import os +import sys +import tkinter as tk +from tkinter import ttk, messagebox, filedialog + +from ccdc.utilities import ApplicationInterface + +from surface_charge_calculator import SurfaceChargeController + + +class SurfaceChargeGUI: + def __init__(self, initial_file_path=None): + self.root = tk.Tk() + self.root.title("Surface Charge Calculator") + try: + photo = tk.PhotoImage(file=os.path.join(os.path.dirname(__file__), 'assets/csd-python-api-logo.png')) + self.root.wm_iconphoto(False, photo) + except FileNotFoundError: + print("Could not find icon file for app.") + except Exception as e: + print("Unable to load icon") + print(e) # This doesn't seem to work with X11 port forwarding 🤷♀️ + # Disable window resizing + self.root.resizable(False, False) + + self.initial_file_path = initial_file_path + self.create_string_file_inputs() + self.create_input_fields() + self.create_buttons() + self.create_treeview() + self.create_directory_selection() + self.configure_grid() # Ensure grid configuration + if self.initial_file_path: + self.handle_initial_file_path(self.initial_file_path) + + def handle_initial_file_path(self, file_path): + """Handles the initial file path by disabling the input fields and setting the file path.""" + self.file_var.set(file_path) # Set the provided file path + self.string_var.set("") # Clear the string input + + # Disable the input fields + self.string_entry.config(state='disabled') + self.file_entry.config(state='readonly') + self.browse_button.config(state='disabled') + + def configure_grid(self): + self.root.grid_rowconfigure(8, weight=1) + self.root.grid_rowconfigure(9, weight=0) + self.root.grid_rowconfigure(10, weight=0) + + self.root.grid_columnconfigure(0, weight=1) + self.root.grid_columnconfigure(1, weight=1) + self.root.grid_columnconfigure(2, weight=1) + self.root.grid_columnconfigure(3, weight=1) + self.root.grid_columnconfigure(4, weight=1) + self.root.grid_columnconfigure(5, weight=1) + self.root.grid_columnconfigure(6, weight=1) + self.root.grid_columnconfigure(7, weight=1) + + def create_string_file_inputs(self): + tk.Label(self.root, text="Structure").grid(row=0, column=0, columnspan=2, sticky='w') + + tk.Label(self.root, text="Refcode:").grid(row=1, column=0, padx=5, pady=5, sticky='e') + self.string_var = tk.StringVar() + self.string_entry = tk.Entry(self.root, textvariable=self.string_var, validate="key", + validatecommand=(self.root.register(self.on_string_input), "%P")) + self.string_entry.grid(row=1, column=1, padx=5, pady=5, columnspan=2, sticky='ew') + + tk.Label(self.root, text="Select File:").grid(row=2, column=0, padx=5, pady=5, sticky='e') + self.file_var = tk.StringVar() + self.file_entry = tk.Entry(self.root, textvariable=self.file_var, state='readonly') + self.file_entry.grid(row=2, column=1, padx=5, pady=5, columnspan=2, sticky='ew') + self.browse_button = tk.Button(self.root, text="Browse", command=self.browse_file) + self.browse_button.grid(row=2, column=3, padx=5, pady=5, sticky='ew') + + def on_string_input(self, input_value): + if input_value.strip(): + self.browse_button.config(state='disabled') + else: + self.browse_button.config(state='normal') + return True + + def create_input_fields(self): + tk.Label(self.root, text="Select hkl and offset").grid(row=3, column=0, columnspan=2, sticky='w') + + input_frame = tk.Frame(self.root) + input_frame.grid(row=4, column=0, columnspan=8, padx=5, pady=5, sticky='ew') + + input_frame.grid_columnconfigure(0, weight=1) + input_frame.grid_columnconfigure(1, weight=1) + input_frame.grid_columnconfigure(2, weight=1) + input_frame.grid_columnconfigure(3, weight=1) + input_frame.grid_columnconfigure(4, weight=1) + input_frame.grid_columnconfigure(5, weight=1) + input_frame.grid_columnconfigure(6, weight=1) + input_frame.grid_columnconfigure(7, weight=1) + + tk.Label(input_frame, text="h:").grid(row=0, column=0, padx=2, pady=5, sticky='e') + tk.Label(input_frame, text="k:").grid(row=0, column=2, padx=2, pady=5, sticky='e') + tk.Label(input_frame, text="l:").grid(row=0, column=4, padx=2, pady=5, sticky='e') + tk.Label(input_frame, text="offset:").grid(row=0, column=6, padx=2, pady=5, sticky='e') + + self.h_var = tk.IntVar() + self.spin_h = tk.Spinbox(input_frame, from_=-9, to=9, width=2, textvariable=self.h_var) + self.spin_h.grid(row=0, column=1, padx=2, pady=5, sticky='ew') + + self.k_var = tk.IntVar() + self.spin_k = tk.Spinbox(input_frame, from_=-9, to=9, width=2, textvariable=self.k_var) + self.spin_k.grid(row=0, column=3, padx=2, pady=5, sticky='ew') + + self.l_var = tk.IntVar() + self.spin_z = tk.Spinbox(input_frame, from_=-9, to=9, width=2, textvariable=self.l_var) + self.spin_z.grid(row=0, column=5, padx=2, pady=5, sticky='ew') + + self.offset_var = tk.DoubleVar() + self.entry_offset = tk.Entry(input_frame, width=10, textvariable=self.offset_var) + self.entry_offset.grid(row=0, column=7, padx=2, pady=5, sticky='ew') + + def create_buttons(self): + self.add_button = tk.Button(self.root, text="Add Surface", command=self.add_combination) + self.add_button.grid(row=5, column=0, columnspan=2, pady=10, sticky='ew') + + self.delete_button = tk.Button(self.root, text="Delete Selected", command=self.delete_combination) + self.delete_button.grid(row=5, column=2, pady=5, sticky='ew') + + self.reset_button = tk.Button(self.root, text="Reset Fields", command=self.reset_fields) + self.reset_button.grid(row=5, column=3, pady=5, sticky='ew') + + self.create_directory_selection() + + def create_directory_selection(self): + tk.Label(self.root, text="Output Directory:").grid(row=9, column=0, padx=5, pady=5, sticky='e') + + self.dir_var = tk.StringVar(value=os.getcwd()) # Default to current working directory + self.dir_entry = tk.Entry(self.root, textvariable=self.dir_var, state='readonly', width=50) + self.dir_entry.grid(row=9, column=1, padx=5, pady=5, columnspan=3, sticky='ew') + + self.browse_dir_button = tk.Button(self.root, text="Browse", command=self.select_directory) + self.browse_dir_button.grid(row=9, column=4, padx=5, pady=5, sticky='ew') + + self.calculate_button = tk.Button(self.root, text="Calculate", command=self.calculate) + self.calculate_button.grid(row=10, column=0, columnspan=5, pady=10, sticky='ew') + + def select_directory(self): + selected_dir = filedialog.askdirectory(initialdir=self.dir_var.get(), title="Select Output Directory") + if selected_dir: + self.dir_var.set(selected_dir) + + def create_treeview(self): + + tk.Label(self.root, text="Current Selections").grid(row=7, column=0, padx=5, pady=5, columnspan=8, + sticky='w') + self.combination_tree = ttk.Treeview(self.root, columns=("h", "k", "l", "Offset"), show='headings') + self.combination_tree.grid(row=8, column=0, columnspan=8, padx=10, pady=10, sticky='nsew') + + self.combination_tree.heading("h", text="h") + self.combination_tree.heading("k", text="k") + self.combination_tree.heading("l", text="l") + self.combination_tree.heading("Offset", text="Offset") + + self.combination_tree.column("h", width=50, anchor=tk.CENTER) + self.combination_tree.column("k", width=50, anchor=tk.CENTER) + self.combination_tree.column("l", width=50, anchor=tk.CENTER) + self.combination_tree.column("Offset", width=100, anchor=tk.CENTER) + + def browse_file(self): + file_path = filedialog.askopenfilename(filetypes=[("mol2 files", "*.mol2")]) + if file_path: + self.file_var.set(file_path) + + def add_combination(self): + try: + h = self.h_var.get() + k = self.k_var.get() + l = self.l_var.get() + if (h, k, l) == (0, 0, 0): + messagebox.showerror("Invalid input", "Please enter valid integers for h, k, l and a float for offset.") + return + offset = self.offset_var.get() + combination = (h, k, l, offset) + if not self.is_duplicate(combination): + self.combination_tree.insert('', tk.END, values=combination) + else: + messagebox.showwarning("Duplicate Entry", "This hkl and offset already exists.") + except tk.TclError: + messagebox.showerror("Invalid input", "Please enter valid integers for h, k, l and a float for offset.") + + def is_duplicate(self, combination): + combination_converted = tuple((str(i) for i in combination)) + for row_id in self.combination_tree.get_children(): + row_values = self.combination_tree.item(row_id, 'values') + if tuple(row_values) == combination_converted: + return True + return False + + def delete_combination(self): + selected_item = self.combination_tree.selection() + if selected_item: + self.combination_tree.delete(selected_item) + else: + messagebox.showwarning("No selection", "Please select a surface to delete.") + + def reset_fields(self): + self.h_var.set(0) + self.k_var.set(0) + self.l_var.set(0) + self.offset_var.set(0.0) + self.string_var.set("") + self.file_var.set("") + self.browse_button.config(state='normal') + + def calculate(self): + string_input = self.string_var.get().strip() + file_input = self.file_var.get().strip() + if not (string_input or file_input): + tk.messagebox.showerror("Input Error", "Please provide a refcode or select a file.") + return + + if not self.combination_tree.get_children(): + tk.messagebox.showerror("Selection Error", "There must be at least one surface in the list.") + return + + items = self.combination_tree.get_children() + data = [] + for item in items: + values = self.combination_tree.item(item, 'values') + try: + h = int(values[0]) + k = int(values[1]) + l = int(values[2]) + offset = float(values[3]) + data.append((h, k, l, offset)) + except ValueError as e: + print(f"Error converting data: {e}") + continue + if string_input: + input_string = string_input # Use string input if available + elif file_input: + input_string = file_input + + output_dir = self.dir_var.get() + + surface_charge_controller = SurfaceChargeController(structure=input_string, output_directory=output_dir, + hkl_and_offsets=data) + surface_charge_controller.calculate_surface_charge() + surface_charge_controller.make_report() + self.root.destroy() + + +if __name__ == "__main__": + if len(sys.argv) > 3 and sys.argv[3].endswith(".m2a"): + mercury = ApplicationInterface() + run_from_mercury = True + input_structure = mercury.input_mol2_file + app = SurfaceChargeGUI(initial_file_path=input_structure) + app.root.mainloop() + else: + app = SurfaceChargeGUI() + app.root.mainloop() diff --git a/scripts/surface_charge/surface_charge_calculator.py b/scripts/surface_charge/surface_charge_calculator.py new file mode 100644 index 0000000..498c891 --- /dev/null +++ b/scripts/surface_charge/surface_charge_calculator.py @@ -0,0 +1,145 @@ +# +# 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 = """ + +
+hkl | +Offset | +Projected Area | +Surface Charge* | +Surface Charge per Projected Area | +
---|---|---|---|---|
{hkl} | +{offset:.2f} | +{self.projected_area[i]:.3f} | +{self.surface_charges[i]:.3f} | +{self.surface_charges_per_area[i]:.4f} | +
*-Surface charge is based on gasteiger partial charges 10.1016/S0040-4039(01)94977-9
+ + + """ + return html