From b7e267202a4f3b434732325193b5714fdbc2853b Mon Sep 17 00:00:00 2001 From: Jack Clegg Date: Fri, 10 Jul 2026 12:13:23 +1000 Subject: [PATCH 1/4] new cif analysis features, minor graphing updates, fix orthonormal in reconstruction --- cx_asap/cxasap.py | 460 ++++++++++++++- .../modules/ADP_analysis.py | 50 +- .../modules/cif_analysis.py | 209 +++++++ .../modules/cif_geometry.py | 544 ++++++++++++++++++ .../modules/cif_read.py | 206 ++++++- .../modules/rotation_planes.py | 31 +- .../modules/structural_analysis.py | 85 +-- .../pipelines/cif_analysis_pipeline.py | 380 ++++++++++++ .../pipelines/variable_position_analysis.py | 82 +-- cx_asap/system_files/crystal_math.py | 91 +++ cx_asap/system_files/parameter.yaml | 45 ++ cx_asap/system_files/utils.py | 88 ++- .../tools/modules/molecule_reconstruction.py | 48 +- 13 files changed, 2063 insertions(+), 256 deletions(-) create mode 100644 cx_asap/post_refinement_analysis/modules/cif_analysis.py create mode 100644 cx_asap/post_refinement_analysis/modules/cif_geometry.py create mode 100644 cx_asap/post_refinement_analysis/pipelines/cif_analysis_pipeline.py diff --git a/cx_asap/cxasap.py b/cx_asap/cxasap.py index c396bb5..bd16479 100644 --- a/cx_asap/cxasap.py +++ b/cx_asap/cxasap.py @@ -113,7 +113,7 @@ import shutil from typing import Union, Tuple -from system_files.utils import Generate, File_Sorter +from system_files.utils import Generate, File_Sorter, format_yaml_error_message from system_files.test_installation import Test from data_reduction.modules.xprep_intensity_compare import Intensity_Compare from data_reduction.modules.XDS_cell_transformation import XDS_Cell_Transformation @@ -134,8 +134,10 @@ from post_refinement_analysis.modules.structural_analysis import Structural_Analysis from post_refinement_analysis.modules.ADP_analysis import ADP_analysis from post_refinement_analysis.modules.centroids import Centroids +from post_refinement_analysis.modules.cif_analysis import CIF_Analysis from post_refinement_analysis.pipelines.rotation_pipeline import Rotation_Pipeline from post_refinement_analysis.pipelines.centroids_pipeline import Centroids_Pipeline +from post_refinement_analysis.pipelines.cif_analysis_pipeline import CIF_Analysis_Pipeline from post_refinement_analysis.pipelines.variable_cif_parameter import ( Variable_Analysis_Pipeline, ) @@ -289,6 +291,9 @@ def yaml_extraction(heading: str) -> dict: "point_geometry_angles", "point_geometry_torsions", "point_geometry_plane_distances", + "rotation_reference_plane", + "rotation_plane_definitions", + "mean_plane_definitions", ] structure_params = [ @@ -365,6 +370,30 @@ def yaml_extraction(heading: str) -> dict: "plane_symmetry": 0, } ] + elif item == "rotation_reference_plane": + yaml_dict[item] = [0, 0, 0] + elif item == "rotation_plane_definitions": + yaml_dict[item] = [ + { + "label": 0, + "plane_atoms": 0, + "plane_symmetry": 0, + } + ] + elif item == "mean_plane_definitions": + yaml_dict[item] = [ + { + "label": 0, + "plane_atoms": 0, + "plane_symmetry": 0, + } + ] + elif item == "cif_input_mode": + yaml_dict[item] = "nested" + elif item == "precombine_cifs": + yaml_dict[item] = False + elif item == "varying_cif_parameter": + yaml_dict[item] = "_diffrn_ambient_temperature" elif item == "reference_plane" or item == "starting_coordinates": yaml_dict[item] = [0, 0, 0] elif item in list_params: @@ -438,10 +467,18 @@ def configuration_check(heading: str) -> Tuple[bool, dict]: with open(yaml_path, "r") as f: try: cfg = yaml.load(f, yaml.FullLoader) - except: - click.echo("Failed to set up config file. Try reconfiguring") + except yaml.YAMLError as error: + click.echo(format_yaml_error_message(yaml_path, error)) exit() + # If varying parameter is left blank/zero, fall back to Data_Block. + if "varying_cif_parameter" in cfg: + varying_value = cfg.get("varying_cif_parameter") + if varying_value in [0, None]: + cfg["varying_cif_parameter"] = "Data_Block" + elif isinstance(varying_value, str) and varying_value.strip() in ["", "0"]: + cfg["varying_cif_parameter"] = "Data_Block" + flag = True for item in yaml_dict.keys(): @@ -3904,6 +3941,419 @@ def pipeline_point_geometry(dependencies, files, configure, run): click.echo("Please select an option. To view options, add --help") +#####------ Module CIF Analysis -----####### + + +@click.command( + "module-cif-analysis", + short_help="analyse CIF folder with optional point/rotation/ADP/structural outputs", +) +@click.option("--dependencies", is_flag=True, help="view the software dependencies") +@click.option("--files", is_flag=True, help="view the required input files") +@click.option("--configure", is_flag=True, help="generate your conf.yaml file") +@click.option("--run", is_flag=True, help="run the code!") +def module_cif_analysis(dependencies, files, configure, run): + """Runs CIF-based analysis for one folder, with optional point-geometry + and rotation-plane calculations directly from CIF atom coordinates. + """ + if dependencies: + click.echo("\nYou do not require any additional software in your path!\n") + elif files: + click.echo("\nYou require the below files:") + click.echo(" - one folder containing one or more .cif files") + click.echo( + " - optional: one .lst file only if you want legacy MPLA fallback rotation analysis" + ) + click.echo("\nThis folder can be located anywhere ") + elif configure: + click.echo("\nWriting a file called conf.yaml in the cx_asap folder...\n") + click.echo("You will need to fill out the parameters.") + click.echo("Descriptions are listed below:") + click.echo( + " - folder_containing_cifs: full path to the folder containing your .cif files" + ) + click.echo( + " - cif_parameters: cif parameters to extract (defaults are usually enough)" + ) + click.echo( + " - atoms_for_analysis: atom labels for structural-analysis filtering" + ) + click.echo( + " - varying_cif_parameter: cif heading used as x-axis, e.g. _diffrn_ambient_temperature" + ) + click.echo( + " if left blank or set to 0, CX-ASAP falls back to Data_Block" + ) + click.echo( + " examples: _diffrn_ambient_temperature, _diffrn_ambient_pressure" + ) + click.echo( + " - reference_unit_cell: optional path to reference .ins for cell-deformation analysis" + ) + click.echo( + " - structural_analysis_bonds/angles/torsions/hbonds: true or false" + ) + click.echo(" - ADP_analysis: true or false") + click.echo( + " - point_geometry_distances: list of distance definitions with keys label, point_1_atoms, point_2_atoms and optional point_n_symmetry" + ) + click.echo( + " - point_geometry_angles: list of angle definitions with keys label, point_1_atoms, point_2_atoms, point_3_atoms and optional point_n_symmetry" + ) + click.echo( + " - point_geometry_torsions: list of torsion definitions with keys label, point_1_atoms..point_4_atoms and optional point_n_symmetry" + ) + click.echo( + " - point_geometry_plane_distances: list of point-plane definitions with keys label, point_atoms, plane_atoms and optional point_symmetry/plane_symmetry" + ) + click.echo( + " - rotation_reference_plane: optional [h,k,l] reference plane for CIF-native rotation angles" + ) + click.echo( + " - rotation_plane_definitions: optional list of planes defined by plane_atoms and optional plane_symmetry" + ) + click.echo( + " - mean_plane_definitions: optional list of planes used for mean-plane/interplane analysis (plane_atoms and optional plane_symmetry)" + ) + click.echo( + " - calculate_interplane_angle: set true to calculate angle between first two mean_plane_definitions planes" + ) + click.echo( + " - mercury_output: set true for 3 dp rounded-centroid companion outputs for point geometry" + ) + click.echo( + " - lst_file_location: optional explicit .lst path (otherwise first .lst in folder is used)" + ) + + fields = yaml_extraction("module-cif-analysis") + yaml_creation(fields) + + elif run: + click.echo("\nChecking to see if experiment configured....\n") + + check, cfg = configuration_check("module-cif-analysis") + + if check == False: + click.echo("Make sure you fill in the configuration file!") + click.echo( + "If you last ran a different code, make sure you reconfigure for the new script!" + ) + click.echo("Re-run configuration for description of each parameter\n") + else: + click.echo("READY TO RUN SCRIPT!\n") + reset_logs() + + distance_defs = cfg.get("point_geometry_distances") or [] + angle_defs = cfg.get("point_geometry_angles") or [] + torsion_defs = cfg.get("point_geometry_torsions") or [] + plane_defs = cfg.get("point_geometry_plane_distances") or [] + + valid_distance_defs = [ + item + for item in distance_defs + if _point_geometry_definition_complete( + item, ["point_1_atoms", "point_2_atoms"] + ) + ] + valid_angle_defs = [ + item + for item in angle_defs + if _point_geometry_definition_complete( + item, ["point_1_atoms", "point_2_atoms", "point_3_atoms"] + ) + ] + valid_torsion_defs = [ + item + for item in torsion_defs + if _point_geometry_definition_complete( + item, + [ + "point_1_atoms", + "point_2_atoms", + "point_3_atoms", + "point_4_atoms", + ], + ) + ] + valid_plane_defs = [ + item + for item in plane_defs + if _point_geometry_definition_complete( + item, ["point_atoms", "plane_atoms"] + ) + ] + rotation_plane_defs = cfg.get("rotation_plane_definitions") or [] + valid_rotation_plane_defs = [ + item + for item in rotation_plane_defs + if _point_geometry_definition_complete(item, ["plane_atoms"]) + ] + mean_plane_defs = cfg.get("mean_plane_definitions") or [] + valid_mean_plane_defs = [ + item + for item in mean_plane_defs + if _point_geometry_definition_complete(item, ["plane_atoms"]) + ] + + if ( + cfg.get("calculate_interplane_angle", False) + and len(valid_mean_plane_defs) == 0 + and len(valid_rotation_plane_defs) > 0 + ): + click.echo( + "No valid mean_plane_definitions found; reusing rotation_plane_definitions for interplane calculation." + ) + valid_mean_plane_defs = valid_rotation_plane_defs + + results_dir = pathlib.Path(cfg["folder_containing_cifs"]) + analysis = CIF_Analysis() + analysis.run( + cfg["folder_containing_cifs"], + str(results_dir), + cfg["cif_parameters"], + cfg["atoms_for_analysis"], + cfg["varying_cif_parameter"], + reference_unit_cell=cfg.get("reference_unit_cell", ""), + structural_analysis_bonds=cfg["structural_analysis_bonds"], + structural_analysis_angles=cfg["structural_analysis_angles"], + structural_analysis_torsions=cfg["structural_analysis_torsions"], + structural_analysis_hbonds=cfg["structural_analysis_hbonds"], + ADP_analysis_enabled=cfg["ADP_analysis"], + point_geometry_distances=valid_distance_defs, + point_geometry_angles=valid_angle_defs, + point_geometry_torsions=valid_torsion_defs, + point_geometry_plane_distances=valid_plane_defs, + mercury_output=cfg.get("mercury_output", False), + reference_plane=cfg.get("rotation_reference_plane", None), + rotation_plane_definitions=valid_rotation_plane_defs, + mean_plane_definitions=valid_mean_plane_defs, + calculate_interplane_angle=cfg.get("calculate_interplane_angle", False), + lst_file_location=cfg.get("lst_file_location", ""), + ) + + copy_logs(str(results_dir)) + + output_message() + + else: + click.echo("Please select an option. To view options, add --help") + + +#####------ Pipeline CIF Analysis -----####### + + +@click.command( + "pipeline-cif-analysis", + short_help="batch CIF analysis with optional point/rotation/ADP/structural outputs", +) +@click.option("--dependencies", is_flag=True, help="view the software dependencies") +@click.option("--files", is_flag=True, help="view the required input files") +@click.option("--configure", is_flag=True, help="generate your conf.yaml file") +@click.option("--run", is_flag=True, help="run the code!") +def pipeline_cif_analysis(dependencies, files, configure, run): + """Runs CIF-based analysis over either nested dataset folders or one flat CIF folder. + Optional point-geometry and rotation-plane calculations run directly from CIF data. + """ + if dependencies: + click.echo("\nYou do not require any additional software in your path!\n") + elif files: + click.echo("\nYou require the below files:") + click.echo(" - experiment_location containing either:") + click.echo(" a) nested dataset folders with .cif files, or") + click.echo(" b) a single flat folder of .cif files") + click.echo( + " - optional: .lst files only if you want legacy MPLA fallback rotation analysis" + ) + click.echo("\nThis folder can be located anywhere ") + elif configure: + click.echo("\nWriting a file called conf.yaml in the cx_asap folder...\n") + click.echo("You will need to fill out the parameters.") + click.echo("Descriptions are listed below:") + click.echo(" - experiment_location: parent folder for CIF analysis") + click.echo( + " - cif_input_mode: 'nested' for subfolders or 'flat' for one CIF folder" + ) + click.echo( + " root-level CIFs are preferred; if both root and nested CIFs exist, root CIFs are used" + ) + click.echo( + " - precombine_cifs: set true to merge discovered CIFs into one combined input file before analysis" + ) + click.echo(" default: false") + click.echo( + " set true when you want one merged input before extraction/analysis" + ) + click.echo( + " precombine outputs are written per run as combined_input.cif and combined_input_sources.txt in the CIF_Analysis/ folder" + ) + click.echo( + " - cif_parameters: cif parameters to extract (defaults are usually enough)" + ) + click.echo( + " - atoms_for_analysis: atom labels for structural-analysis filtering" + ) + click.echo( + " - varying_cif_parameter: cif heading used as x-axis, e.g. _diffrn_ambient_temperature" + ) + click.echo( + " if left blank or set to 0, CX-ASAP falls back to Data_Block" + ) + click.echo( + " examples: _diffrn_ambient_temperature, _diffrn_ambient_pressure" + ) + click.echo( + " - reference_unit_cell: optional path to reference .ins for cell-deformation analysis" + ) + click.echo( + " - structural_analysis_bonds/angles/torsions/hbonds: true or false" + ) + click.echo(" - ADP_analysis: true or false") + click.echo( + " - point_geometry_distances: list of distance definitions with keys label, point_1_atoms, point_2_atoms and optional point_n_symmetry" + ) + click.echo( + " - point_geometry_angles: list of angle definitions with keys label, point_1_atoms, point_2_atoms, point_3_atoms and optional point_n_symmetry" + ) + click.echo( + " - point_geometry_torsions: list of torsion definitions with keys label, point_1_atoms..point_4_atoms and optional point_n_symmetry" + ) + click.echo( + " - point_geometry_plane_distances: list of point-plane definitions with keys label, point_atoms, plane_atoms and optional point_symmetry/plane_symmetry" + ) + click.echo( + " - rotation_reference_plane: optional [h,k,l] reference plane for CIF-native rotation angles" + ) + click.echo( + " - rotation_plane_definitions: optional list of planes defined by plane_atoms and optional plane_symmetry" + ) + click.echo( + " - mean_plane_definitions: optional list of planes used for mean-plane/interplane analysis (plane_atoms and optional plane_symmetry)" + ) + click.echo( + " - calculate_interplane_angle: set true to calculate angle between first two mean_plane_definitions planes" + ) + click.echo( + " - mercury_output: set true for 3 dp rounded-centroid companion outputs for point geometry" + ) + + fields = yaml_extraction("pipeline-cif-analysis") + yaml_creation(fields) + + elif run: + click.echo("\nChecking to see if experiment configured....\n") + + check, cfg = configuration_check("pipeline-cif-analysis") + + if check == False: + click.echo("Make sure you fill in the configuration file!") + click.echo( + "If you last ran a different code, make sure you reconfigure for the new script!" + ) + click.echo("Re-run configuration for description of each parameter\n") + else: + click.echo("READY TO RUN SCRIPT!\n") + reset_logs() + + distance_defs = cfg.get("point_geometry_distances") or [] + angle_defs = cfg.get("point_geometry_angles") or [] + torsion_defs = cfg.get("point_geometry_torsions") or [] + plane_defs = cfg.get("point_geometry_plane_distances") or [] + + valid_distance_defs = [ + item + for item in distance_defs + if _point_geometry_definition_complete( + item, ["point_1_atoms", "point_2_atoms"] + ) + ] + valid_angle_defs = [ + item + for item in angle_defs + if _point_geometry_definition_complete( + item, ["point_1_atoms", "point_2_atoms", "point_3_atoms"] + ) + ] + valid_torsion_defs = [ + item + for item in torsion_defs + if _point_geometry_definition_complete( + item, + [ + "point_1_atoms", + "point_2_atoms", + "point_3_atoms", + "point_4_atoms", + ], + ) + ] + valid_plane_defs = [ + item + for item in plane_defs + if _point_geometry_definition_complete( + item, ["point_atoms", "plane_atoms"] + ) + ] + rotation_plane_defs = cfg.get("rotation_plane_definitions") or [] + valid_rotation_plane_defs = [ + item + for item in rotation_plane_defs + if _point_geometry_definition_complete(item, ["plane_atoms"]) + ] + mean_plane_defs = cfg.get("mean_plane_definitions") or [] + valid_mean_plane_defs = [ + item + for item in mean_plane_defs + if _point_geometry_definition_complete(item, ["plane_atoms"]) + ] + + if ( + cfg.get("calculate_interplane_angle", False) + and len(valid_mean_plane_defs) == 0 + and len(valid_rotation_plane_defs) > 0 + ): + click.echo( + "No valid mean_plane_definitions found; reusing rotation_plane_definitions for interplane calculation." + ) + valid_mean_plane_defs = valid_rotation_plane_defs + + pipe = CIF_Analysis_Pipeline() + results_dir = pipe.create_numbered_results_directory( + cfg["experiment_location"], "CIF_Analysis" + ) + + pipe.run( + cfg["experiment_location"], + str(results_dir), + cfg["cif_input_mode"], + cfg["cif_parameters"], + cfg["atoms_for_analysis"], + cfg["varying_cif_parameter"], + reference_unit_cell=cfg.get("reference_unit_cell", ""), + structural_analysis_bonds=cfg["structural_analysis_bonds"], + structural_analysis_angles=cfg["structural_analysis_angles"], + structural_analysis_torsions=cfg["structural_analysis_torsions"], + structural_analysis_hbonds=cfg["structural_analysis_hbonds"], + ADP_analysis_enabled=cfg["ADP_analysis"], + point_geometry_distances=valid_distance_defs, + point_geometry_angles=valid_angle_defs, + point_geometry_torsions=valid_torsion_defs, + point_geometry_plane_distances=valid_plane_defs, + mercury_output=cfg.get("mercury_output", False), + reference_plane=cfg.get("rotation_reference_plane", None), + rotation_plane_definitions=valid_rotation_plane_defs, + mean_plane_definitions=valid_mean_plane_defs, + calculate_interplane_angle=cfg.get("calculate_interplane_angle", False), + precombine_cifs=cfg.get("precombine_cifs", False), + ) + + copy_logs(str(results_dir)) + + output_message() + + else: + click.echo("Please select an option. To view options, add --help") + + #######------Pipeline varying parameter ---------###### """This pipeline will analyse .cif files for a dynamic experiment where one @@ -5057,6 +5507,8 @@ def pipeline_AS_Brute_individual(dependencies, files, configure, run): pipeline_shelxt_auto, module_point_geometry, pipeline_point_geometry, + module_cif_analysis, + pipeline_cif_analysis, ] if BadOS == True: @@ -5107,6 +5559,8 @@ def pipeline_AS_Brute_individual(dependencies, files, configure, run): cli.add_command(pipeline_rotation_planes) cli.add_command(module_point_geometry) cli.add_command(pipeline_point_geometry) + cli.add_command(module_cif_analysis) + cli.add_command(pipeline_cif_analysis) cli.add_command(pipeline_position_analysis) cli.add_command(pipeline_AS_Brute) cli.add_command(module_molecule_reconstruction) diff --git a/cx_asap/post_refinement_analysis/modules/ADP_analysis.py b/cx_asap/post_refinement_analysis/modules/ADP_analysis.py index 92895cf..771a173 100755 --- a/cx_asap/post_refinement_analysis/modules/ADP_analysis.py +++ b/cx_asap/post_refinement_analysis/modules/ADP_analysis.py @@ -211,48 +211,14 @@ def analyse_data(self, csv_file: str, cell_data: str) -> None: values_unsorted, vectors = np.linalg.eig(BG) - test = values_unsorted.tolist() - - values = sorted(test) - - values.reverse() - - positions = {} - - for i_1, i in enumerate(values_unsorted): - for i_2, j in enumerate(values): - if i == j: - positions[i_1] = i_2 - - vector_1 = np.array( - [ - [ - vectors[0][positions[0]], - vectors[1][positions[0]], - vectors[2][positions[0]], - ] - ] - ) - - vector_2 = np.array( - [ - [ - vectors[0][positions[1]], - vectors[1][positions[1]], - vectors[2][positions[1]], - ] - ] - ) - - vector_3 = np.array( - [ - [ - vectors[0][positions[2]], - vectors[1][positions[2]], - vectors[2][positions[2]], - ] - ] - ) + # Numerical noise can produce tiny imaginary components. + # Sort by real part and coerce near-real values back to float. + order = np.argsort(np.real(values_unsorted))[::-1] + values = [float(np.real(values_unsorted[i])) for i in order] + + vector_1 = np.array([np.real_if_close(vectors[:, order[0]]).astype(float)]) + vector_2 = np.array([np.real_if_close(vectors[:, order[1]]).astype(float)]) + vector_3 = np.array([np.real_if_close(vectors[:, order[2]]).astype(float)]) try: principle_A.append(math.sqrt(values[0] / (2 * (math.pi**2)))) diff --git a/cx_asap/post_refinement_analysis/modules/cif_analysis.py b/cx_asap/post_refinement_analysis/modules/cif_analysis.py new file mode 100644 index 0000000..e924870 --- /dev/null +++ b/cx_asap/post_refinement_analysis/modules/cif_analysis.py @@ -0,0 +1,209 @@ +#!/usr/bin/env python3 + +################################################################################################### +# --------------------------------------CX-ASAP: cif_analysis--------------------------------------# +# ---Authors: Amy J. Thompson, Kate M. Smith, Daniel J. Eriksson, Jack K. Clegg & Jason R. Price---# +# -----------------------------------Python Implementation by AJT----------------------------------# +# -----------------------------------Project Design by JRP and JKC---------------------------------# +# --------------------------------Valuable Coding Support by KMS & DJE-----------------------------# +################################################################################################### + +from system_files.utils import Config +from post_refinement_analysis.modules.cif_read import CIF_Read +from post_refinement_analysis.modules.structural_analysis import Structural_Analysis +from post_refinement_analysis.modules.ADP_analysis import ADP_analysis +from post_refinement_analysis.modules.cell_analysis import Cell_Deformation +from post_refinement_analysis.modules.cif_geometry import CIF_Geometry +from post_refinement_analysis.modules.rotation_planes import Rotation +import pathlib +import os +import logging + + +class CIF_Analysis: + """Coordinates CIF parameter extraction and optional analysis modules.""" + + def __init__(self, test_mode: bool = False) -> None: + """Initialise shared configuration for one analysis run.""" + + self.test_mode = test_mode + + config = Config(self.test_mode) + self.cfg = config.cfg + self.sys = config.sys + self.conf_path = config.conf_path + self.sys_path = config.sys_path + + @staticmethod + def _find_first_file(location: str, suffix: str) -> str: + """Return the first file matching suffix from a file or folder path.""" + + base = pathlib.Path(location) + + if base.is_file(): + if base.suffix.lower() == suffix.lower(): + return str(base) + return "" + + if not base.exists() or not base.is_dir(): + return "" + + for item in base.iterdir(): + if item.is_file() and item.suffix.lower() == suffix.lower(): + return str(item) + return "" + + def run( + self, + cif_location: str, + results_directory: str, + cif_parameters: list, + atoms_for_analysis: list, + varying_cif_parameter: str, + reference_unit_cell: str = "", + structural_analysis_bonds: bool = False, + structural_analysis_angles: bool = False, + structural_analysis_torsions: bool = False, + structural_analysis_hbonds: bool = False, + ADP_analysis_enabled: bool = False, + point_geometry_distances: list = None, + point_geometry_angles: list = None, + point_geometry_torsions: list = None, + point_geometry_plane_distances: list = None, + mercury_output: bool = False, + reference_plane: list = None, + rotation_plane_definitions: list = None, + mean_plane_definitions: list = None, + calculate_interplane_angle: bool = False, + lst_file_location: str = "", + ) -> None: + """Run CIF extraction plus structural/ADP/cell/geometry analysis outputs. + + The method always writes results into ``results_directory``. Geometry and + rotation analyses are optional and are enabled by non-empty definitions. + Legacy LST-based rotation fallback is used only when a reference plane is + provided and no CIF-native rotation plane definitions are configured. + """ + + os.chdir(results_directory) + + point_geometry_distances = point_geometry_distances or [] + point_geometry_angles = point_geometry_angles or [] + point_geometry_torsions = point_geometry_torsions or [] + point_geometry_plane_distances = point_geometry_plane_distances or [] + rotation_plane_definitions = rotation_plane_definitions or [] + mean_plane_definitions = mean_plane_definitions or [] + + cif_data = CIF_Read(self.test_mode) + cif_data.configure(cif_parameters) + cif_data.get_data( + cif_location, + structural_analysis_bonds, + structural_analysis_angles, + structural_analysis_torsions, + structural_analysis_hbonds, + ADP_analysis_enabled, + varying_cif_parameter, + ) + cif_data.data_output() + + geometry = Structural_Analysis(self.test_mode) + bond_file = ( + "Bond_Lengths.csv" + if structural_analysis_bonds and pathlib.Path("Bond_Lengths.csv").exists() + else False + ) + angle_file = ( + "Bond_Angles.csv" + if structural_analysis_angles and pathlib.Path("Bond_Angles.csv").exists() + else False + ) + torsion_file = ( + "Bond_Torsions.csv" + if structural_analysis_torsions + and pathlib.Path("Bond_Torsions.csv").exists() + else False + ) + hbond_file = ( + "HBond_details.csv" + if structural_analysis_hbonds and pathlib.Path("HBond_details.csv").exists() + else False + ) + + geometry.import_and_analyse( + bond_file, + angle_file, + torsion_file, + hbond_file, + atoms_for_analysis, + results_directory, + varying_parameter=varying_cif_parameter, + ) + + if reference_unit_cell not in ["", 0, None]: + cell = Cell_Deformation(self.test_mode) + cell.import_data("CIF_Parameters.csv", reference_unit_cell) + cell.calculate_deformations() + cell.quality_analysis( + varying_cif_parameter, + { + "R1": cell.df["_refine_ls_R_factor_gt"], + "Rint": cell.df["_diffrn_reflns_av_R_equivalents"], + "Completeness": cell.df["_diffrn_measured_fraction_theta_full"], + }, + varying_cif_parameter, + ) + cell.graphical_analysis(varying_cif_parameter, varying_cif_parameter) + + if ADP_analysis_enabled and pathlib.Path("ADPs.csv").exists(): + adp_obj = ADP_analysis(self.test_mode) + adp_obj.analyse_data("ADPs.csv", "CIF_Parameters.csv") + + run_cif_geometry = ( + len(point_geometry_distances) > 0 + or len(point_geometry_angles) > 0 + or len(point_geometry_torsions) > 0 + or len(point_geometry_plane_distances) > 0 + or ( + reference_plane not in [None, 0, "", [0], [0, 0, 0]] + and len(rotation_plane_definitions) > 0 + ) + or (calculate_interplane_angle and len(mean_plane_definitions) > 0) + ) + + if run_cif_geometry: + cif_geometry = CIF_Geometry(self.test_mode) + cif_geometry.run( + cif_location, + results_directory, + varying_cif_parameter, + point_geometry_distances, + point_geometry_angles, + point_geometry_torsions, + point_geometry_plane_distances, + mercury_output=mercury_output, + reference_plane=reference_plane, + rotation_plane_definitions=rotation_plane_definitions, + mean_plane_definitions=mean_plane_definitions, + calculate_interplane_angle=calculate_interplane_angle, + ) + + if lst_file_location in ["", 0, None]: + lst_file_location = self._find_first_file(cif_location, ".lst") + + # Legacy fallback: retain old MPLA-on-LST behaviour when no CIF plane definitions are supplied. + if ( + reference_plane not in [None, 0, "", [0], [0, 0, 0]] + and len(rotation_plane_definitions) == 0 + ): + if lst_file_location not in ["", 0, None]: + rot = Rotation(self.test_mode) + rot.configure(reference_plane) + rot.analysis(lst_file_location, 1, results_directory) + if calculate_interplane_angle: + rot.analyse_interplane_angle(lst_file_location, 1, results_directory) + else: + logging.warning( + __name__ + + " : Rotation-plane analysis requested but no CIF rotation_plane_definitions or .lst file was found" + ) diff --git a/cx_asap/post_refinement_analysis/modules/cif_geometry.py b/cx_asap/post_refinement_analysis/modules/cif_geometry.py new file mode 100644 index 0000000..8ba4865 --- /dev/null +++ b/cx_asap/post_refinement_analysis/modules/cif_geometry.py @@ -0,0 +1,544 @@ +#!/usr/bin/env python3 + +################################################################################################### +# --------------------------------------CX-ASAP: cif_geometry--------------------------------------# +# ---Authors: Amy J. Thompson, Kate M. Smith, Daniel J. Eriksson, Jack K. Clegg & Jason R. Price---# +# -----------------------------------Python Implementation by AJT----------------------------------# +# -----------------------------------Project Design by JRP and JKC---------------------------------# +# --------------------------------Valuable Coding Support by KMS & DJE-----------------------------# +################################################################################################### + +from CifFile import ReadCif +from post_refinement_analysis.modules.centroids import Centroids +from system_files.crystal_math import ( + fractional_to_cartesian, + best_fit_plane_normal, + cartesian_plane_normal_to_fractional, + angle_between_vectors, +) +import pathlib +import pandas as pd +import logging +import numpy as np + + +class CIF_Geometry: + """Performs symmetry-aware point and plane analysis directly from CIF files.""" + + def __init__(self, test_mode: bool = False) -> None: + """Initialise CIF geometry engine and centroid helper.""" + + self.test_mode = test_mode + self.point_engine = Centroids(self.test_mode) + + @staticmethod + def _parse_float(raw) -> float: + """Parses a CIF numeric value, handling values with uncertainties.""" + + if isinstance(raw, (int, float)): + return float(raw) + + raw_str = str(raw).strip() + if raw_str in ["", ".", "?"]: + raise ValueError("Missing numeric CIF value") + + if "(" in raw_str: + raw_str = raw_str.split("(")[0] + + return float(raw_str) + + @staticmethod + def _read_cif_files(cif_location: str) -> list: + """Discover CIF inputs from a file path or root-first folder layout.""" + + base = pathlib.Path(cif_location) + + if base.is_file(): + if base.suffix.lower() == ".cif": + return [base.resolve()] + return [] + + if not base.exists() or not base.is_dir(): + return [] + + def _is_results_folder(folder_name: str) -> bool: + name = folder_name.strip().lower() + blocked_exact = { + "cif_analysis", + "geometry_analysis", + "refinement_statistics", + "results", + "analysis", + "ref", + "failed_autoprocessing", + } + if name in blocked_exact: + return True + if name.startswith("_"): + return True + return False + + root_files = [item for item in sorted(base.glob("*.cif")) if item.is_file()] + + files = [] + + # One folder level down only + for child in sorted(base.iterdir()): + if not child.is_dir(): + continue + if _is_results_folder(child.name): + continue + files.extend( + [item for item in sorted(child.glob("*.cif")) if item.is_file()] + ) + + if len(root_files) > 0: + if len(files) > 0: + logging.warning( + __name__ + + " : Mixed CIF layout detected (root-level and nested CIFs). " + + "Using root-level CIFs only; nested CIFs will be ignored." + ) + return [item.resolve() for item in root_files] + + # Deduplicate while preserving deterministic order + seen = set() + unique_files = [] + for item in files: + key = str(item.resolve()) + if key in seen: + continue + seen.add(key) + unique_files.append(item) + + return unique_files + + def _extract_structure_data(self, block) -> dict: + """Extract unit-cell and fractional atom coordinates from one CIF block.""" + + cell_keys = [ + "_cell_length_a", + "_cell_length_b", + "_cell_length_c", + "_cell_angle_alpha", + "_cell_angle_beta", + "_cell_angle_gamma", + ] + cell = [self._parse_float(block[item]) for item in cell_keys] + + labels = block["_atom_site_label"] + x_vals = block["_atom_site_fract_x"] + y_vals = block["_atom_site_fract_y"] + z_vals = block["_atom_site_fract_z"] + + if not isinstance(labels, list): + labels = [labels] + x_vals = [x_vals] + y_vals = [y_vals] + z_vals = [z_vals] + + coords = {} + for label, x_raw, y_raw, z_raw in zip(labels, x_vals, y_vals, z_vals): + try: + coords[str(label).upper()] = [ + self._parse_float(x_raw), + self._parse_float(y_raw), + self._parse_float(z_raw), + ] + except ValueError: + continue + + return { + "block": block, + "coords": coords, + "cell": cell, + } + + @staticmethod + def _append_rows(csv_name: str, rows: list, results_directory: str) -> None: + """Write collected row dictionaries to a CSV when rows are present.""" + + if len(rows) == 0: + return + + df = pd.DataFrame(rows) + output_path = pathlib.Path(results_directory) / csv_name + df.to_csv(output_path, index=None) + + @staticmethod + def _definition_has_centroid(definition: dict, keys: list, engine: Centroids) -> bool: + """Check whether any definition key uses centroid-style atom syntax.""" + + for key in keys: + if engine._point_uses_centroid(definition.get(key, [])): + return True + return False + + @staticmethod + def _valid_reference_plane(reference_plane: list) -> bool: + """Validate a reference plane is a numeric [h, k, l] list.""" + + if not isinstance(reference_plane, list) or len(reference_plane) != 3: + return False + try: + [float(item) for item in reference_plane] + except (TypeError, ValueError): + return False + return True + + def _rotation_plane_angle( + self, + coords: dict, + cell_params: list, + definition: dict, + reference_plane: list, + ) -> tuple: + """Compute acute angle between a fitted plane normal and reference plane.""" + + normal_frac = self._plane_normal_fractional(coords, cell_params, definition) + ref_frac = np.array(reference_plane, dtype=float) + angle = angle_between_vectors(normal_frac, ref_frac, fold_to_acute=True) + + return angle, normal_frac + + def _plane_normal_fractional( + self, + coords: dict, + cell_params: list, + definition: dict, + ) -> np.ndarray: + """Fit a plane from atoms and return its normal in fractional coordinates.""" + + atoms = definition.get("plane_atoms", []) + symmetry = definition.get("plane_symmetry") or None + frac_points = self.point_engine._extract_positions(coords, atoms, symmetry) + + if len(frac_points) < 3: + raise ValueError("Need at least 3 atoms to define rotation plane") + + cart_points = np.array( + [fractional_to_cartesian(item, cell_params) for item in frac_points], + dtype=float, + ) + normal_cart = best_fit_plane_normal(cart_points) + normal_frac = cartesian_plane_normal_to_fractional(normal_cart, cell_params) + return normal_frac + + def run( + self, + cif_location: str, + results_directory: str, + varying_parameter: str, + point_geometry_distances: list = None, + point_geometry_angles: list = None, + point_geometry_torsions: list = None, + point_geometry_plane_distances: list = None, + mercury_output: bool = False, + reference_plane: list = None, + rotation_plane_definitions: list = None, + mean_plane_definitions: list = None, + calculate_interplane_angle: bool = False, + ) -> None: + """Run CIF-native point geometry and rotation-plane analyses. + + For each CIF block, this writes CSV outputs for any configured distance, + angle, torsion, point-plane, rotation-plane, and optional interplane + calculations. Mercury-style rounded-centroid outputs are generated when + ``mercury_output`` is enabled and centroid definitions are present. + """ + + point_geometry_distances = point_geometry_distances or [] + point_geometry_angles = point_geometry_angles or [] + point_geometry_torsions = point_geometry_torsions or [] + point_geometry_plane_distances = point_geometry_plane_distances or [] + rotation_plane_definitions = rotation_plane_definitions or [] + mean_plane_definitions = mean_plane_definitions or [] + + files = self._read_cif_files(cif_location) + if len(files) == 0: + logging.warning(__name__ + " : No .cif files found for CIF geometry analysis") + return + + rows_distance = [] + rows_angle = [] + rows_torsion = [] + rows_plane = [] + rows_distance_mercury = [] + rows_angle_mercury = [] + rows_torsion_mercury = [] + rows_plane_mercury = [] + rows_rotation = [] + rows_interplane = [] + + structure_number = 0 + for cif_file in files: + cif_obj = ReadCif(str(cif_file)) + for block_name in cif_obj.keys(): + structure_number += 1 + + try: + structure = self._extract_structure_data(cif_obj[block_name]) + except Exception as error: + logging.warning( + __name__ + + " : Failed to extract structure from " + + str(cif_file) + + " block " + + str(block_name) + + " due to " + + str(error) + ) + continue + + block = structure["block"] + coords = structure["coords"] + cell = structure["cell"] + self.point_engine.cell_params = cell + + base_row = { + "Structure": structure_number, + "CIF_File": cif_file.stem, + "Data_Block": block_name, + } + try: + base_row[varying_parameter] = self._parse_float(block[varying_parameter]) + except Exception: + pass + + distance_row = dict(base_row) + angle_row = dict(base_row) + torsion_row = dict(base_row) + plane_row = dict(base_row) + distance_row_mercury = dict(base_row) + angle_row_mercury = dict(base_row) + torsion_row_mercury = dict(base_row) + plane_row_mercury = dict(base_row) + + for index, definition in enumerate(point_geometry_distances): + if not isinstance(definition, dict): + continue + label = definition.get("label", f"Distance_{index + 1}") + distance_row[label] = self.point_engine.point_distance( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + ) + if mercury_output and self._definition_has_centroid( + definition, + ["point_1_atoms", "point_2_atoms"], + self.point_engine, + ): + distance_row_mercury[label] = self.point_engine.point_distance( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + centroid_round_dp=3, + ) + + for index, definition in enumerate(point_geometry_angles): + if not isinstance(definition, dict): + continue + label = definition.get("label", f"Angle_{index + 1}") + angle_row[label] = self.point_engine.point_angle( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_3_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + definition.get("point_3_symmetry") or None, + ) + if mercury_output and self._definition_has_centroid( + definition, + ["point_1_atoms", "point_2_atoms", "point_3_atoms"], + self.point_engine, + ): + angle_row_mercury[label] = self.point_engine.point_angle( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_3_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + definition.get("point_3_symmetry") or None, + centroid_round_dp=3, + ) + + for index, definition in enumerate(point_geometry_torsions): + if not isinstance(definition, dict): + continue + label = definition.get("label", f"Torsion_{index + 1}") + torsion_row[label] = self.point_engine.point_torsion( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_3_atoms", []), + definition.get("point_4_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + definition.get("point_3_symmetry") or None, + definition.get("point_4_symmetry") or None, + ) + if mercury_output and self._definition_has_centroid( + definition, + [ + "point_1_atoms", + "point_2_atoms", + "point_3_atoms", + "point_4_atoms", + ], + self.point_engine, + ): + torsion_row_mercury[label] = self.point_engine.point_torsion( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_3_atoms", []), + definition.get("point_4_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + definition.get("point_3_symmetry") or None, + definition.get("point_4_symmetry") or None, + centroid_round_dp=3, + ) + + for index, definition in enumerate(point_geometry_plane_distances): + if not isinstance(definition, dict): + continue + label = definition.get("label", f"Point_Plane_Distance_{index + 1}") + plane_row[label] = self.point_engine.point_plane_distance( + coords, + definition.get("point_atoms", []), + definition.get("plane_atoms", []), + definition.get("point_symmetry") or None, + definition.get("plane_symmetry") or None, + ) + if mercury_output and self._definition_has_centroid( + definition, + ["point_atoms"], + self.point_engine, + ): + plane_row_mercury[label] = self.point_engine.point_plane_distance( + coords, + definition.get("point_atoms", []), + definition.get("plane_atoms", []), + definition.get("point_symmetry") or None, + definition.get("plane_symmetry") or None, + centroid_round_dp=3, + ) + + if len(distance_row.keys()) > len(base_row.keys()): + rows_distance.append(distance_row) + if len(angle_row.keys()) > len(base_row.keys()): + rows_angle.append(angle_row) + if len(torsion_row.keys()) > len(base_row.keys()): + rows_torsion.append(torsion_row) + if len(plane_row.keys()) > len(base_row.keys()): + rows_plane.append(plane_row) + + if mercury_output: + if len(distance_row_mercury.keys()) > len(base_row.keys()): + rows_distance_mercury.append(distance_row_mercury) + if len(angle_row_mercury.keys()) > len(base_row.keys()): + rows_angle_mercury.append(angle_row_mercury) + if len(torsion_row_mercury.keys()) > len(base_row.keys()): + rows_torsion_mercury.append(torsion_row_mercury) + if len(plane_row_mercury.keys()) > len(base_row.keys()): + rows_plane_mercury.append(plane_row_mercury) + + if self._valid_reference_plane(reference_plane): + rotation_row = dict(base_row) + for index, definition in enumerate(rotation_plane_definitions): + if not isinstance(definition, dict): + continue + + label = definition.get("label", f"MPLA_{index + 1}_Rotation_Angle") + try: + angle, normal_frac = self._rotation_plane_angle( + coords, + cell, + definition, + reference_plane, + ) + except ValueError as error: + logging.warning(__name__ + " : " + str(error)) + continue + + rotation_row[label] = angle + + if len(rotation_row.keys()) > len(base_row.keys()): + rows_rotation.append(rotation_row) + + if calculate_interplane_angle: + mean_normals = [] + mean_labels = [] + for index, definition in enumerate(mean_plane_definitions): + if not isinstance(definition, dict): + continue + + try: + mean_normals.append( + self._plane_normal_fractional(coords, cell, definition) + ) + mean_labels.append( + definition.get("label", f"Mean_Plane_{index + 1}") + ) + except ValueError as error: + logging.warning(__name__ + " : " + str(error)) + continue + + if len(mean_normals) >= 2: + interplane = angle_between_vectors( + mean_normals[0], mean_normals[1], fold_to_acute=True + ) + rows_interplane.append( + { + "Structure": structure_number, + "CIF_File": cif_file.stem, + "Data_Block": block_name, + "Mean Plane 1": mean_labels[0], + "Mean Plane 2": mean_labels[1], + "Interplane Angle": interplane, + } + ) + elif len(mean_plane_definitions) > 0: + logging.warning( + __name__ + + " : calculate_interplane_angle requires at least two valid mean planes" + ) + + self._append_rows("point_geometry_distances.csv", rows_distance, results_directory) + self._append_rows("point_geometry_angles.csv", rows_angle, results_directory) + self._append_rows("point_geometry_torsions.csv", rows_torsion, results_directory) + self._append_rows( + "point_geometry_plane_distances.csv", rows_plane, results_directory + ) + + if mercury_output: + self._append_rows( + "point_geometry_distances_mercury.csv", + rows_distance_mercury, + results_directory, + ) + self._append_rows( + "point_geometry_angles_mercury.csv", + rows_angle_mercury, + results_directory, + ) + self._append_rows( + "point_geometry_torsions_mercury.csv", + rows_torsion_mercury, + results_directory, + ) + self._append_rows( + "point_geometry_plane_distances_mercury.csv", + rows_plane_mercury, + results_directory, + ) + + self._append_rows("rotation_angles.csv", rows_rotation, results_directory) + self._append_rows("interplane_angles.csv", rows_interplane, results_directory) diff --git a/cx_asap/post_refinement_analysis/modules/cif_read.py b/cx_asap/post_refinement_analysis/modules/cif_read.py index 4c05a48..8f77eb4 100755 --- a/cx_asap/post_refinement_analysis/modules/cif_read.py +++ b/cx_asap/post_refinement_analysis/modules/cif_read.py @@ -22,6 +22,8 @@ class CIF_Read: + """Reads CIF inputs and extracts tabular data for downstream analysis.""" + def __init__(self, test_mode: bool = False) -> None: """Initialises the class @@ -54,6 +56,7 @@ def __init__(self, test_mode: bool = False) -> None: self.successful_positions = [] self.results = {} self.errors = {} + self._cif_cache = {} # Sets these to 0 to reset from previous runs @@ -81,6 +84,7 @@ def configure(self, search_items: list) -> None: # Pulls parameters from the configuration file as necessary, and uses it to set up an empty dataframe self.search_items = search_items + self._cif_cache = {} for item in self.search_items: self.results[item] = [] @@ -93,6 +97,83 @@ def configure(self, search_items: list) -> None: self.temp_df = pd.DataFrame() self.adp_data = pd.DataFrame() + def _read_cif(self, cif_file: pathlib.Path): + """Reads a CIF once per file path during a run and reuses it.""" + + cache_key = str(pathlib.Path(cif_file).resolve()) + if cache_key not in self._cif_cache: + self._cif_cache[cache_key] = ReadCif(str(pathlib.Path(cif_file))) + return self._cif_cache[cache_key] + + @staticmethod + def _read_cif_files(location: str) -> list: + """Collects CIF files from a file path, root, or one level below. + + Root-first policy: + - If root contains CIFs, use only those. + - If root has no CIFs, scan one level down while skipping likely + CX-ASAP results folders. + """ + + base = pathlib.Path(location) + + if base.is_file(): + if base.suffix.lower() == ".cif": + return [base.resolve()] + return [] + + if not base.exists() or not base.is_dir(): + return [] + + def _is_results_folder(folder_name: str) -> bool: + name = folder_name.strip().lower() + blocked_exact = { + "cif_analysis", + "geometry_analysis", + "refinement_statistics", + "results", + "analysis", + "ref", + "failed_autoprocessing", + } + if name in blocked_exact: + return True + if name.startswith("_"): + return True + return False + + root_files = [item for item in sorted(base.glob("*.cif")) if item.is_file()] + + files = [] + for child in sorted(base.iterdir()): + if not child.is_dir(): + continue + if _is_results_folder(child.name): + continue + files.extend( + [item for item in sorted(child.glob("*.cif")) if item.is_file()] + ) + + if len(root_files) > 0: + if len(files) > 0: + logging.warning( + __name__ + + " : Mixed CIF layout detected (root-level and nested CIFs). " + + "Using root-level CIFs only; nested CIFs will be ignored." + ) + return [item.resolve() for item in root_files] + + seen = set() + unique_files = [] + for item in files: + key = str(item.resolve()) + if key in seen: + continue + seen.add(key) + unique_files.append(item) + + return unique_files + def parameter_tidy(self, raw: str, item: str) -> None: """This tidies the output from the CIF and separates @@ -181,14 +262,15 @@ def get_data( # This function searches through all of the folders in the current working directory for a cif file - self.tree_browse = Directory_Browse(pathlib.Path(location), self.test_mode) - - self.tree_browse.enter_directory_multiple(pathlib.Path(location), ".cif") + cif_files = self._read_cif_files(location) + if len(cif_files) == 0: + logging.warning(__name__ + " : No .cif files found in " + str(location)) + return # For all found cif_files: - for index, item in enumerate(self.tree_browse.item_files): - cif_file = self.tree_browse.item_files[index].absolute() + for cif_file in cif_files: + cif_obj = self._read_cif(cif_file) # extracts the desired cif parameters, as well as how many structures per cif and which positions were successful @@ -196,12 +278,20 @@ def get_data( temp_data, structures_in_cif_tmp, successful_positions_tmp, - ) = self.data_harvest(cif_file, self.search_items, varying_parameter) + ) = self.data_harvest( + cif_file, self.search_items, varying_parameter, cif_obj=cif_obj + ) self.structural_analysis( - cif_file, bonds, angles, torsions, hbonds, varying_parameter + cif_file, + bonds, + angles, + torsions, + hbonds, + varying_parameter, + cif_obj=cif_obj, ) - self.adp_analysis(cif_file, adp) + self.adp_analysis(cif_file, adp, varying_parameter, cif_obj=cif_obj) # self.data = self.data.append(temp_data) self.data = pd.concat([self.data, temp_data]) @@ -229,6 +319,7 @@ def structural_analysis( torsions: bool = False, hbonds: bool = False, varying_parameter: str = "_diffrn_ambient_temperature", + cif_obj=None, ) -> None: """Extracts structural information from CIF @@ -281,7 +372,9 @@ def structural_analysis( temp_data_bonds, structures_in_cif_tmp_bonds, successful_positions_tmp_bonds, - ) = self.data_harvest(cif_file, bond_paras, varying_parameter) + ) = self.data_harvest( + cif_file, bond_paras, varying_parameter, cif_obj=cif_obj + ) # self.bond_data = self.bond_data.append(temp_data_bonds) self.bond_data = pd.concat([self.bond_data, temp_data_bonds]) if angles == True: @@ -289,7 +382,9 @@ def structural_analysis( temp_data_angles, structures_in_cif_tmp_angles, successful_positions_tmp_angles, - ) = self.data_harvest(cif_file, angle_paras, varying_parameter) + ) = self.data_harvest( + cif_file, angle_paras, varying_parameter, cif_obj=cif_obj + ) # self.angle_data = self.angle_data.append(temp_data_angles) self.angle_data = pd.concat([self.angle_data, temp_data_angles]) if torsions == True: @@ -297,7 +392,9 @@ def structural_analysis( temp_data_torsions, structures_in_cif_tmp_torsions, successful_positions_tmp_torsions, - ) = self.data_harvest(cif_file, torsion_paras, varying_parameter) + ) = self.data_harvest( + cif_file, torsion_paras, varying_parameter, cif_obj=cif_obj + ) # self.torsion_data = self.torsion_data.append(temp_data_torsions) self.torsion_data = pd.concat([self.torsion_data, temp_data_torsions]) if hbonds == True: @@ -305,7 +402,9 @@ def structural_analysis( temp_data_hbonds, structures_in_cif_tmp_hbonds, successful_positions_tmp_hbonds, - ) = self.data_harvest(cif_file, hbond_paras, varying_parameter) + ) = self.data_harvest( + cif_file, hbond_paras, varying_parameter, cif_obj=cif_obj + ) # self.hbond_data = self.hbond_data.append(temp_data_hbonds) self.hbond_data = pd.concat([self.hbond_data, temp_data_hbonds]) @@ -314,6 +413,7 @@ def adp_analysis( cif_file: str, adp: bool = False, varying_parameter: str = "_diffrn_ambient_temperature", + cif_obj=None, ) -> None: """Extracts ADP information from CIF @@ -337,7 +437,9 @@ def adp_analysis( temp_data_adps, structures_in_cif_tmp_adps, successful_positions_tmp_adps, - ) = self.data_harvest(cif_file, adps, varying_parameter) + ) = self.data_harvest( + cif_file, adps, varying_parameter, cif_obj=cif_obj + ) # self.adp_data = self.adp_data.append(temp_data_adps) self.adp_data = pd.concat([self.adp_data, temp_data_adps]) @@ -346,6 +448,7 @@ def data_harvest( cif_file: str, search_items: list, varying_parameter: str = "_diffrn_ambient_temperature", + cif_obj=None, ) -> Tuple["pd.DataFrame", int, list]: """Extracts all other desired parameters from CIF @@ -373,7 +476,7 @@ def data_harvest( # Use of the PyCifRW library for easy parsing of CIF Files - cif = ReadCif(cif_file.name) + cif = cif_obj if cif_obj is not None else self._read_cif(pathlib.Path(cif_file)) # Identifies datablocks within the CIF File @@ -396,9 +499,14 @@ def data_harvest( self.cif_list = [] structure_analysis_counter[item] = [] - for experiment in self.data_blocks: + for block_index, experiment in enumerate(self.data_blocks): try: - raw = cif[experiment][item] + if item == "Data_Block": + raw = str(experiment) + elif item == "Structure": + raw = block_index + 1 + else: + raw = cif[experiment][item] except: logging.critical("Failed to find " + item + " in " + cif_file.stem) print("Critical Failure - see error log for details") @@ -442,21 +550,65 @@ def data_harvest( self.results[item] += [numbers_to_multiply[index]] * i self.errors[item] += [errors_to_multiply[index]] * i - self.temp_df["CIF_File"] = self.cif_list + # Build a robust per-datablock repeat profile from all harvested parameters. + counters = list(structure_analysis_counter.values()) + target_counter = [] + if len(counters) != 0: + n_blocks = len(counters[0]) + for block_index in range(n_blocks): + target_counter.append( + max(counter_list[block_index] for counter_list in counters) + ) + target_len = len(self.cif_list) + if len(target_counter) != 0: + target_len = max(target_len, sum(target_counter)) for para in search_items: - if len(self.results[para]) != len(self.cif_list): - self.temp_df = pd.DataFrame() + target_len = max(target_len, len(self.results[para])) - for para in search_items: - self.temp_df[para] = self.results[para] - self.temp_df[para + "_error"] = self.errors[para] + def _expand_or_pad(values: list, target: int, counts: list) -> list: + if len(values) == target: + return values - if len(self.temp_df) != len(self.cif_list): - ( - self.temp_df["CIF_File"], - self.temp_df["Data_Block"], - ) = self.generate_cif_list(self.temp_df, test_val) + if len(values) == len(counts): + expanded = [] + for idx, count in enumerate(counts): + expanded += [values[idx]] * count + if len(expanded) == target: + return expanded + + if len(values) == 1 and target > 1: + return values * target + + if len(values) < target: + return values + [None] * (target - len(values)) + + return values[:target] + + self.temp_df = pd.DataFrame(index=range(target_len)) + + if target_len == len(self.cif_list): + self.temp_df["CIF_File"] = self.cif_list + self.temp_df["Data_Block"] = list(self.data_blocks) + else: + expanded_cif = [] + expanded_blocks = [] + block_names = list(self.data_blocks) + if len(target_counter) != 0: + for index, count in enumerate(target_counter): + expanded_cif += [self.cif_list[index]] * count + expanded_blocks += [block_names[index]] * count + + expanded_cif = _expand_or_pad(expanded_cif, target_len, target_counter) + expanded_blocks = _expand_or_pad(expanded_blocks, target_len, target_counter) + self.temp_df["CIF_File"] = expanded_cif + self.temp_df["Data_Block"] = expanded_blocks + + for para in search_items: + values = _expand_or_pad(self.results[para], target_len, target_counter) + errors = _expand_or_pad(self.errors[para], target_len, target_counter) + self.temp_df[para] = values + self.temp_df[para + "_error"] = errors return self.temp_df, number_of_structures, self.data_blocks diff --git a/cx_asap/post_refinement_analysis/modules/rotation_planes.py b/cx_asap/post_refinement_analysis/modules/rotation_planes.py index e114ed1..8253c14 100755 --- a/cx_asap/post_refinement_analysis/modules/rotation_planes.py +++ b/cx_asap/post_refinement_analysis/modules/rotation_planes.py @@ -11,7 +11,11 @@ # ----------Required Modules----------# from system_files.utils import Nice_YAML_Dumper, Config -from system_files.crystal_math import orthonorm_matrix +from system_files.crystal_math import ( + cartesian_plane_normal_to_fractional, + angle_between_vectors, + reciprocal_orthonorm_matrix, +) from post_refinement_analysis.modules.lst_read import LST_Read import pathlib import os @@ -133,30 +137,13 @@ def calculate_planes(self, data: list, ref_plane: list, ref_values: list) -> lis logging.warning(__name__ + " : No MPLA planes found in .lst file") return [] - M = orthonorm_matrix(self.ref_values) - M_star = np.linalg.inv(M) - - # Convert reference plane vector to fractional space - ref = np.array( - [[self.ref_plane[0], self.ref_plane[1], self.ref_plane[2]]], dtype=float - ) - ref_frac = np.dot(ref, M_star) + M_star = reciprocal_orthonorm_matrix(self.ref_values) + ref_frac = np.dot(np.array(self.ref_plane, dtype=float), M_star) angles = [] for normal in plane_normals: - cart = np.array([[normal[0], normal[1], normal[2]]], dtype=float) - frac = np.dot(cart, M_star) - - angle = float( - np.degrees( - np.arccos( - np.dot(frac, ref_frac.T) - / (np.linalg.norm(frac) * np.linalg.norm(ref_frac)) - ) - )[0][0] - ) - if 180 - angle < 90: - angle = 180 - angle + frac = cartesian_plane_normal_to_fractional(normal, self.ref_values) + angle = angle_between_vectors(frac, ref_frac, fold_to_acute=True) angles.append(angle) return angles diff --git a/cx_asap/post_refinement_analysis/modules/structural_analysis.py b/cx_asap/post_refinement_analysis/modules/structural_analysis.py index 871bae2..f706be0 100755 --- a/cx_asap/post_refinement_analysis/modules/structural_analysis.py +++ b/cx_asap/post_refinement_analysis/modules/structural_analysis.py @@ -261,27 +261,9 @@ def structural_analysis( important_df = pd.DataFrame() if df.empty == False: - # Separates out important atoms by looking for them in any column and merging into one dataframe - - for item in atoms_for_analysis: - temp_df = df[df.eq(item).any(axis=1)] - important_df = pd.concat([important_df, temp_df], axis=0) - - # Need to make a new column of the indices - the above code will give you double ups if both atoms in a bond are "important" - - # When concatanating the dataframes, it keeps the indicies of the original dataframes - - # Meaning... doubleup indices AND out of order indices - - # Both are bad, so make a new column of the indicies, use it to drop duplicates, and then resets the index of the important dataframe - - important_df["index_2"] = list(important_df.index) - - important_df = important_df.drop_duplicates(subset=["index_2"]) - - important_df = important_df.drop(["index_2"], axis=1) - - important_df = important_df.reset_index(drop=True) + # One-pass mask is equivalent to the previous per-atom concat/deduplicate flow. + important_mask = df.isin(atoms_for_analysis).any(axis=1) + important_df = df[important_mask].copy().reset_index(drop=True) # Making symmetry equivalent bonds (ie same atom 1 and atom 2) distinguishable @@ -316,39 +298,12 @@ def structural_analysis( logging.info("Something went weird.") dup = important_df.duplicated(["Joined", varying_parameter], keep=False) - - # The below function counts the number of each group of duplicates - - list_dup = important_df.pivot_table( - columns=["Joined", varying_parameter], aggfunc="size" - ).to_dict() - - counter = 0 - - new_column = [] - - # This appends a suffix to each duplicated bond in the joined column based on how many there are - - for j, i in enumerate(dup): - bond = important_df["Joined"][j] - - # NOTE HERE VARIABLE 'TEMPERATURE' CAN BE ANY PARAMETER BUT I DIDN'T WANT TO CHANGE WHOLE CODE - - temperature = important_df[varying_parameter][j] - - if i == True: - temp = important_df["Joined"][j] - new_column.append(important_df["Joined"][j] + "_" + str(counter)) - counter += 1 - else: - new_column.append(important_df["Joined"][j]) - try: - if counter == list_dup[(bond, temperature)]: - counter = 0 - except UnboundLocalError: - pass - - important_df["Joined"] = new_column + dup_order = important_df.groupby(["Joined", varying_parameter]).cumcount() + important_df.loc[dup, "Joined"] = ( + important_df.loc[dup, "Joined"] + + "_" + + dup_order[dup].astype(str) + ) important_df.to_csv(prefix + "_Important_" + file_name, index=None) @@ -384,7 +339,7 @@ def structural_analysis( os.chdir(folder_name) for item in discrete_atoms: - separated_df = df[df.eq(item).any(axis=1)] + separated_df = df[df["Joined"] == item] separated_df.to_csv( structure_type + "_" + str(item) + ".csv", index=None ) @@ -418,19 +373,17 @@ def structural_analysis( if flexible == True: x_data.append(list(g["number"])) x_unit = "Structure number" - x_data.append(list(g[varying_parameter])) - if column_names[-4] == "_diffrn_ambient_temperature": - x_unit = "Temperature (K)" - elif ( - flexible == False - and column_names[-4] != "_diffrn_ambient_temperature" - ): - x_unit = varying_parameter + else: + x_data.append(list(g[varying_parameter])) + if column_names[-4] == "_diffrn_ambient_temperature": + x_unit = "Temperature (K)" + else: + x_unit = varying_parameter if structure_type == "Hbonds": try: graph.single_scatter_graph( - x_data[0], + x_data, da_data, x_unit, r"D$\cdots$A Distance ($\AA$)", @@ -445,7 +398,7 @@ def structural_analysis( try: graph.single_scatter_graph( - x_data[0], + x_data, angle_data, x_unit, "D-H$\\cdots$A Angle ($^\\circ$)", @@ -460,7 +413,7 @@ def structural_analysis( else: try: graph.single_scatter_graph( - x_data[0], + x_data, y_data, x_unit, y_unit, diff --git a/cx_asap/post_refinement_analysis/pipelines/cif_analysis_pipeline.py b/cx_asap/post_refinement_analysis/pipelines/cif_analysis_pipeline.py new file mode 100644 index 0000000..4350360 --- /dev/null +++ b/cx_asap/post_refinement_analysis/pipelines/cif_analysis_pipeline.py @@ -0,0 +1,380 @@ +#!/usr/bin/env python3 + +################################################################################################### +# ----------------------------------CX-ASAP: cif_analysis_pipeline---------------------------------# +# ---Authors: Amy J. Thompson, Kate M. Smith, Daniel J. Eriksson, Jack K. Clegg & Jason R. Price---# +# -----------------------------------Python Implementation by AJT----------------------------------# +# -----------------------------------Project Design by JRP and JKC---------------------------------# +# --------------------------------Valuable Coding Support by KMS & DJE-----------------------------# +################################################################################################### + +from system_files.utils import Config +from system_files.utils import Grapher +from post_refinement_analysis.modules.cif_analysis import CIF_Analysis +from CifFile import ReadCif +import pathlib +import logging +import pandas as pd + + +class CIF_Analysis_Pipeline: + """Pipeline wrapper for running CIF analysis on directories or files.""" + + def __init__(self, test_mode: bool = False) -> None: + """Initialise pipeline configuration and runtime options.""" + + self.test_mode = test_mode + + config = Config(self.test_mode) + self.cfg = config.cfg + self.sys = config.sys + self.conf_path = config.conf_path + self.sys_path = config.sys_path + + def create_numbered_results_directory( + self, base_directory: str, folder_name: str + ) -> pathlib.Path: + """Create the next numbered results folder under ``base_directory``.""" + + base_path = pathlib.Path(base_directory) + results_root = base_path / folder_name + results_root.mkdir(exist_ok=True) + + existing_numbers = [] + for item in results_root.iterdir(): + if item.is_dir(): + try: + existing_numbers.append(int(item.name)) + except ValueError: + continue + + next_number = max(existing_numbers) + 1 if existing_numbers else 1 + results_path = results_root / str(next_number) + results_path.mkdir() + return results_path + + def _dataset_directories(self, experiment_location: str, cif_input_mode: str) -> list: + """Return dataset folders for flat or nested experiment layouts.""" + + mode = str(cif_input_mode).strip().lower() + exp = pathlib.Path(experiment_location) + + if mode == "flat": + return [exp] + + nested_dirs = [item for item in exp.iterdir() if item.is_dir()] + if len(nested_dirs) == 0: + return [exp] + + return nested_dirs + + @staticmethod + def _is_results_folder(folder_name: str) -> bool: + """Identify folders that should be excluded from CIF discovery.""" + + name = folder_name.strip().lower() + blocked_exact = { + "cif_analysis", + "geometry_analysis", + "refinement_statistics", + "results", + "analysis", + "ref", + "failed_autoprocessing", + } + if name in blocked_exact: + return True + if name.startswith("_"): + return True + return False + + def _discover_cif_files(self, location: str) -> list: + """Root-first CIF discovery with one-level fallback.""" + + base = pathlib.Path(location) + + if base.is_file(): + if base.suffix.lower() == ".cif": + return [base.resolve()] + return [] + + if not base.exists() or not base.is_dir(): + return [] + + root_files = [item for item in sorted(base.glob("*.cif")) if item.is_file()] + + nested_files = [] + for child in sorted(base.iterdir()): + if not child.is_dir(): + continue + if self._is_results_folder(child.name): + continue + nested_files.extend( + [item for item in sorted(child.glob("*.cif")) if item.is_file()] + ) + + if len(root_files) > 0: + if len(nested_files) > 0: + logging.warning( + __name__ + + " : Mixed CIF layout detected (root-level and nested CIFs). " + + "Using root-level CIFs only; nested CIFs will be ignored." + ) + return [item.resolve() for item in root_files] + + seen = set() + unique_files = [] + for item in nested_files: + key = str(item.resolve()) + if key in seen: + continue + seen.add(key) + unique_files.append(item) + + return unique_files + + @staticmethod + def _write_combined_cif(cif_files: list, combined_path: pathlib.Path) -> None: + """Merge multiple CIF files into one combined CIF file.""" + + with open(combined_path, "w") as out_file: + for cif_file in cif_files: + cif_obj = ReadCif(str(cif_file)) + out_file.write(cif_obj.WriteOut()) + + @staticmethod + def _write_combined_manifest(cif_files: list, manifest_path: pathlib.Path) -> None: + """Write source-file provenance list for a combined CIF input.""" + + with open(manifest_path, "w") as manifest: + for cif_file in cif_files: + manifest.write(str(cif_file) + "\n") + + @staticmethod + def _combine_dataset_csvs(results_directory: str, csv_name: str) -> "pd.DataFrame | None": + """Combines one CSV from all dataset_* folders into a single pipeline-level CSV.""" + + root = pathlib.Path(results_directory) + root_csv = root / csv_name + + # If analysis already wrote directly to the root folder, reuse that file. + if root_csv.exists(): + try: + return pd.read_csv(root_csv) + except Exception as error: + logging.warning( + __name__ + + " : Could not read " + + str(root_csv) + + " due to " + + str(error) + ) + + frames = [] + + for dataset_dir in sorted(root.glob("dataset_*")): + csv_path = dataset_dir / csv_name + if not csv_path.exists(): + continue + + try: + df = pd.read_csv(csv_path) + except Exception as error: + logging.warning( + __name__ + + " : Could not read " + + str(csv_path) + + " due to " + + str(error) + ) + continue + + if len(df) == 0: + continue + + df["Dataset"] = dataset_dir.name + frames.append(df) + + if len(frames) == 0: + return None + + combined = pd.concat(frames, ignore_index=True) + combined.to_csv(root / csv_name, index=None) + return combined + + @staticmethod + def _plot_combined_csv( + df: "pd.DataFrame | None", + varying_cif_parameter: str, + y_axis_title: str, + graph_title: str, + figure_name: str, + results_directory: str, + ) -> None: + """Creates a points-style scatter plot for a combined CSV output.""" + + if df is None or len(df) == 0: + return + + if varying_cif_parameter in df.columns: + x_col = varying_cif_parameter + x_title = varying_cif_parameter + elif "Structure" in df.columns: + x_col = "Structure" + x_title = "Structure Number" + else: + return + + metadata_cols = {"Structure", "CIF_File", "Data_Block", "Dataset", x_col} + y_cols = [col for col in df.columns if col not in metadata_cols] + y_cols = [col for col in y_cols if pd.api.types.is_numeric_dtype(df[col])] + + if not y_cols: + return + + x_data = list(df[x_col]) + y_data = [list(df[col]) for col in y_cols] + + graph = Grapher() + graph.single_scatter_graph( + x_data, + y_data, + x_title, + y_axis_title, + graph_title, + str(pathlib.Path(results_directory) / figure_name), + y_series_title=y_cols if len(y_cols) > 1 else None, + ) + + def run( + self, + experiment_location: str, + results_directory: str, + cif_input_mode: str, + cif_parameters: list, + atoms_for_analysis: list, + varying_cif_parameter: str, + reference_unit_cell: str = "", + structural_analysis_bonds: bool = False, + structural_analysis_angles: bool = False, + structural_analysis_torsions: bool = False, + structural_analysis_hbonds: bool = False, + ADP_analysis_enabled: bool = False, + point_geometry_distances: list = None, + point_geometry_angles: list = None, + point_geometry_torsions: list = None, + point_geometry_plane_distances: list = None, + mercury_output: bool = False, + reference_plane: list = None, + rotation_plane_definitions: list = None, + mean_plane_definitions: list = None, + calculate_interplane_angle: bool = False, + precombine_cifs: bool = False, + ) -> None: + """Execute a CIF analysis pipeline run for a chosen experiment location. + + When ``precombine_cifs`` is true and multiple inputs are discovered, a + combined CIF and source manifest are created in the run output folder, + then analysed as a single input. + """ + + analyser = CIF_Analysis(self.test_mode) + + logging.info( + __name__ + + " : Running CIF analysis pipeline with mode=" + + str(cif_input_mode) + + " at root location=" + + str(experiment_location) + ) + + analysis_location = str(experiment_location) + + if precombine_cifs: + cif_files = self._discover_cif_files(str(experiment_location)) + if len(cif_files) == 0: + logging.warning( + __name__ + + " : precombine_cifs is true but no CIF files were discovered at " + + str(experiment_location) + ) + elif len(cif_files) == 1: + analysis_location = str(cif_files[0]) + logging.info( + __name__ + + " : precombine_cifs requested with one CIF; analysing file directly: " + + analysis_location + ) + else: + combined_path = pathlib.Path(results_directory) / "combined_input.cif" + manifest_path = pathlib.Path(results_directory) / "combined_input_sources.txt" + self._write_combined_cif(cif_files, combined_path) + self._write_combined_manifest(cif_files, manifest_path) + analysis_location = str(combined_path) + logging.info( + __name__ + + " : precombine_cifs created combined input: " + + str(combined_path) + ) + + # Run once at the selected analysis location to keep output layout + # aligned with regular CIF analysis (single results folder). + analyser.run( + analysis_location, + str(results_directory), + cif_parameters, + atoms_for_analysis, + varying_cif_parameter, + reference_unit_cell=reference_unit_cell, + structural_analysis_bonds=structural_analysis_bonds, + structural_analysis_angles=structural_analysis_angles, + structural_analysis_torsions=structural_analysis_torsions, + structural_analysis_hbonds=structural_analysis_hbonds, + ADP_analysis_enabled=ADP_analysis_enabled, + point_geometry_distances=point_geometry_distances, + point_geometry_angles=point_geometry_angles, + point_geometry_torsions=point_geometry_torsions, + point_geometry_plane_distances=point_geometry_plane_distances, + mercury_output=mercury_output, + reference_plane=reference_plane, + rotation_plane_definitions=rotation_plane_definitions, + mean_plane_definitions=mean_plane_definitions, + calculate_interplane_angle=calculate_interplane_angle, + ) + + # Consolidate angle outputs at pipeline root to match points-pipeline UX. + point_angles_df = self._combine_dataset_csvs( + results_directory, "point_geometry_angles.csv" + ) + self._plot_combined_csv( + point_angles_df, + varying_cif_parameter, + "Angle($^\\circ$)", + "Point Geometry Angles", + "point_geometry_angles.png", + results_directory, + ) + + rotation_angles_df = self._combine_dataset_csvs( + results_directory, "rotation_angles.csv" + ) + self._plot_combined_csv( + rotation_angles_df, + varying_cif_parameter, + "Angle($^\\circ$)", + "Rotation Angles", + "rotation_angles.png", + results_directory, + ) + + interplane_angles_df = self._combine_dataset_csvs( + results_directory, "interplane_angles.csv" + ) + self._plot_combined_csv( + interplane_angles_df, + varying_cif_parameter, + "Angle($^\\circ$)", + "Interplane Angles", + "interplane_angles.png", + results_directory, + ) diff --git a/cx_asap/post_refinement_analysis/pipelines/variable_position_analysis.py b/cx_asap/post_refinement_analysis/pipelines/variable_position_analysis.py index 1aee431..948ee29 100755 --- a/cx_asap/post_refinement_analysis/pipelines/variable_position_analysis.py +++ b/cx_asap/post_refinement_analysis/pipelines/variable_position_analysis.py @@ -189,7 +189,7 @@ def analyse_data( "Rint": item["_diffrn_reflns_av_R_equivalents"], "Completeness": item["_diffrn_measured_fraction_theta_full"], }, - "Distance($\mu$m)", + "Distance($\\mu$m)", df=item, ) @@ -197,7 +197,7 @@ def analyse_data( cell.graphical_analysis( "Distance", - "Distance($\mu$m)", + "Distance($\\mu$m)", "Cell Deformation - " + discrete_cif_names[index], "cell_parameters_" + discrete_cif_names[index] + ".png", "axis_deformation_" + discrete_cif_names[index] + ".png", @@ -209,19 +209,14 @@ def analyse_data( if adps != False: adp_df = pd.read_csv("ADPs.csv") - discrete_cif_names_bond = list(dict.fromkeys(bond_df["CIF_File"])) - separated_by_cif_bond = [] - - for item in discrete_cif_names_bond: - condition = bond_df["CIF_File"] == item - separated_by_cif_bond.append(bond_df[condition]) - - for index, item in enumerate(separated_by_cif_bond): - item.to_csv("ADPs_" + discrete_cif_names_bond[index] + ".csv") + for index, (cif_name, item) in enumerate( + adp_df.groupby("CIF_File", sort=False) + ): + item.to_csv("ADPs_" + cif_name + ".csv", index=None) adp_object = ADP_analysis(self.test_mode) adp_object.analyse_data( - "ADPs_" + discrete_cif_names_bond[index] + ".csv", + "ADPs_" + cif_name + ".csv", discrete_cif_names[index] + "_parameters.csv", ) @@ -231,18 +226,13 @@ def analyse_data( if bonds != False: bond_df = pd.read_csv("Bond_Lengths.csv") - discrete_cif_names_bond = list(dict.fromkeys(bond_df["CIF_File"])) - separated_by_cif_bond = [] - - for item in discrete_cif_names_bond: - condition = bond_df["CIF_File"] == item - separated_by_cif_bond.append(bond_df[condition]) - - for index, item in enumerate(separated_by_cif_bond): - item.to_csv("Bond_Lengths_" + discrete_cif_names_bond[index] + ".csv", index=None) + for index, (cif_name, item) in enumerate( + bond_df.groupby("CIF_File", sort=False) + ): + item.to_csv("Bond_Lengths_" + cif_name + ".csv", index=None) geometry.import_and_analyse( - "Bond_Lengths_" + discrete_cif_names_bond[index] + ".csv", + "Bond_Lengths_" + cif_name + ".csv", False, False, False, @@ -254,19 +244,14 @@ def analyse_data( if angles != False: angle_df = pd.read_csv("Bond_Angles.csv") - discrete_cif_names_angle = list(dict.fromkeys(angle_df["CIF_File"])) - separated_by_cif_angle = [] - - for item in discrete_cif_names_angle: - condition = angle_df["CIF_File"] == item - separated_by_cif_angle.append(angle_df[condition]) - - for index, item in enumerate(separated_by_cif_angle): - item.to_csv("Bond_Angles_" + discrete_cif_names_angle[index] + ".csv", index=None) + for index, (cif_name, item) in enumerate( + angle_df.groupby("CIF_File", sort=False) + ): + item.to_csv("Bond_Angles_" + cif_name + ".csv", index=None) geometry.import_and_analyse( False, - "Bond_Angles_" + discrete_cif_names_angle[index] + ".csv", + "Bond_Angles_" + cif_name + ".csv", False, False, atoms_for_analysis, @@ -277,23 +262,15 @@ def analyse_data( if torsions != False: torsion_df = pd.read_csv("Bond_Torsions.csv") - discrete_cif_names_torsion = list(dict.fromkeys(torsion_df["CIF_File"])) - separated_by_cif_torsion = [] - - for item in discrete_cif_names_torsion: - condition = torsion_df["CIF_File"] == item - separated_by_cif_torsion.append(torsion_df[condition]) - - for index, item in enumerate(separated_by_cif_torsion): - item.to_csv( - "Bond_Torsions_" + discrete_cif_names_torsion[index] + ".csv", - index=None, - ) + for index, (cif_name, item) in enumerate( + torsion_df.groupby("CIF_File", sort=False) + ): + item.to_csv("Bond_Torsions_" + cif_name + ".csv", index=None) geometry.import_and_analyse( False, False, - "Bond_Torsions_" + discrete_cif_names_torsion[index] + ".csv", + "Bond_Torsions_" + cif_name + ".csv", False, atoms_for_analysis, location, @@ -302,21 +279,16 @@ def analyse_data( ) if hbonds != False: hbond_df = pd.read_csv("HBond_details.csv") - discrete_cif_names_hbond = list(dict.fromkeys(hbond_df["CIF_File"])) - separated_by_cif_hbond = [] - - for item in discrete_cif_names_hbond: - condition = hbond_df["CIF_File"] == item - separated_by_cif_hbond.append(hbond_df[condition]) - - for index, item in enumerate(separated_by_cif_hbond): - item.to_csv("HBond_details_" + discrete_cif_names_hbond[index] + ".csv", index=None) + for index, (cif_name, item) in enumerate( + hbond_df.groupby("CIF_File", sort=False) + ): + item.to_csv("HBond_details_" + cif_name + ".csv", index=None) geometry.import_and_analyse( False, False, False, - "HBond_details_" + discrete_cif_names_hbond[index] + ".csv", + "HBond_details_" + cif_name + ".csv", atoms_for_analysis, location, str(index + 1), diff --git a/cx_asap/system_files/crystal_math.py b/cx_asap/system_files/crystal_math.py index cfc3ef3..a66ad30 100644 --- a/cx_asap/system_files/crystal_math.py +++ b/cx_asap/system_files/crystal_math.py @@ -245,3 +245,94 @@ def point_to_plane_distance( normal_unit = normal / normal_norm return float(abs(np.dot(p - centroid, normal_unit))) + + +def reciprocal_orthonorm_matrix(cell_params: list) -> np.ndarray: + """Returns the reciprocal-space transform matrix M* = inv(M). + + This is useful for converting vectors between Cartesian and fractional + representations in reciprocal space. + """ + + M = orthonorm_matrix(cell_params) + return np.linalg.inv(M) + + +def cartesian_plane_normal_to_fractional( + normal_cart: "np.ndarray | list", cell_params: list +) -> np.ndarray: + """Converts a Cartesian plane normal into fractional representation. + + Args: + normal_cart: plane normal components in Cartesian basis + cell_params: [a, b, c, alpha, beta, gamma] + + Returns: + np.ndarray: length-3 fractional vector + """ + + n_cart = np.array(normal_cart, dtype=float) + M_star = reciprocal_orthonorm_matrix(cell_params) + return np.dot(n_cart, M_star) + + +def angle_between_vectors( + vector_1: "np.ndarray | list", + vector_2: "np.ndarray | list", + fold_to_acute: bool = False, +) -> float: + """Calculates the angle in degrees between two vectors. + + Args: + vector_1: first vector + vector_2: second vector + fold_to_acute: if True, folds obtuse results into [0, 90] + + Returns: + float: angle in degrees + """ + + v1 = np.array(vector_1, dtype=float) + v2 = np.array(vector_2, dtype=float) + + n1 = np.linalg.norm(v1) + n2 = np.linalg.norm(v2) + if n1 == 0.0 or n2 == 0.0: + raise ValueError("Zero-length vector in angle calculation") + + cos_theta = np.dot(v1, v2) / (n1 * n2) + cos_theta = np.clip(cos_theta, -1.0, 1.0) + angle = float(np.degrees(np.arccos(cos_theta))) + + if fold_to_acute and angle > 90.0: + angle = 180.0 - angle + + return angle + + +def best_fit_plane_normal( + plane_points: "list[np.ndarray] | np.ndarray | list", +) -> np.ndarray: + """Returns a unit normal vector for a best-fit Cartesian plane. + + Args: + plane_points: array-like collection of 3D Cartesian points + + Returns: + np.ndarray: length-3 unit normal vector + """ + + plane = np.array(plane_points, dtype=float) + if len(plane) < 3: + raise ValueError("Need at least 3 points to define a plane") + + centroid = np.mean(plane, axis=0) + centered = plane - centroid + _, _, vh = np.linalg.svd(centered) + normal = vh[-1] + + normal_norm = np.linalg.norm(normal) + if normal_norm == 0.0: + raise ValueError("Could not resolve plane normal") + + return normal / normal_norm diff --git a/cx_asap/system_files/parameter.yaml b/cx_asap/system_files/parameter.yaml index 7611a73..d3e33a2 100644 --- a/cx_asap/system_files/parameter.yaml +++ b/cx_asap/system_files/parameter.yaml @@ -327,6 +327,28 @@ module-point-geometry: - point_geometry_plane_distances - mercury_output +module-cif-analysis: + - folder_containing_cifs + - cif_parameters + - atoms_for_analysis + - varying_cif_parameter + - reference_unit_cell + - structural_analysis_bonds + - structural_analysis_angles + - structural_analysis_torsions + - structural_analysis_hbonds + - ADP_analysis + - point_geometry_distances + - point_geometry_angles + - point_geometry_torsions + - point_geometry_plane_distances + - mercury_output + - rotation_reference_plane + - rotation_plane_definitions + - mean_plane_definitions + - calculate_interplane_angle + - lst_file_location + pipeline-centroids: - experiment_location - centroid_1_atoms @@ -342,6 +364,29 @@ pipeline-point-geometry: - point_geometry_plane_distances - mercury_output +pipeline-cif-analysis: + - experiment_location + - cif_input_mode + - precombine_cifs + - cif_parameters + - atoms_for_analysis + - varying_cif_parameter + - reference_unit_cell + - structural_analysis_bonds + - structural_analysis_angles + - structural_analysis_torsions + - structural_analysis_hbonds + - ADP_analysis + - point_geometry_distances + - point_geometry_angles + - point_geometry_torsions + - point_geometry_plane_distances + - mercury_output + - rotation_reference_plane + - rotation_plane_definitions + - mean_plane_definitions + - calculate_interplane_angle + module-structural-analysis: - atoms_for_analysis - bond_data diff --git a/cx_asap/system_files/utils.py b/cx_asap/system_files/utils.py index 8a3f55c..c33ca82 100755 --- a/cx_asap/system_files/utils.py +++ b/cx_asap/system_files/utils.py @@ -1157,9 +1157,18 @@ def single_scatter_graph( s (list): marker sizes for multiple series """ + x_is_per_series = ( + isinstance(x, list) + and len(x) == len(y) + and len(x) > 0 + and isinstance(x[0], list) + ) + for index, item in enumerate(y): + x_series = x[index] if x_is_per_series else x + if type(item) != float and type(item) != int: - if len(x) != len(item): + if len(x_series) != len(item): logging.info( __name__ + " : Possible error with plotting structural changes. Check the structures in the output CIF for unreasonable structures." @@ -1175,17 +1184,17 @@ def single_scatter_graph( # all the temperatures hadn't been edited yet - to_repeat = x[0] + to_repeat = x_series[0] - x = [to_repeat] * len(item) + x_series = [to_repeat] * len(item) if colour == None and y_series_title != None: - plt.scatter(x, item, label=y_series_title[index]) + plt.scatter(x_series, item, label=y_series_title[index]) elif y_series_title == None: - plt.scatter(x, y) + plt.scatter(x_series, item) else: plt.scatter( - x, + x_series, item, c=colour[index], marker=marker[index], @@ -1195,7 +1204,7 @@ def single_scatter_graph( ) else: - plt.scatter(x, y) + plt.scatter(x_series, item) plt.xlabel(x_title, fontsize=12) plt.ylabel(y_title, fontsize=12) @@ -1847,6 +1856,40 @@ def write_line_break(self, data: dict = None) -> None: super().write_line_break() +def format_yaml_error_message(file_path: pathlib.Path, error: Exception) -> str: + """Formats YAML parsing errors with location and actionable hints.""" + + file_text = str(file_path) + message = [f"Failed to parse YAML file: {file_text}"] + + problem = getattr(error, "problem", None) + if problem: + message.append(f"Problem: {problem}") + + mark = getattr(error, "problem_mark", None) + if mark is not None: + message.append(f"Location: line {mark.line + 1}, column {mark.column + 1}") + + hint = ( + "Hint: check indentation and list formatting, and ensure key/value pairs " + "use a space after ':'." + ) + + problem_text = str(problem).lower() if problem else "" + if "mapping values are not allowed here" in problem_text: + hint = ( + "Hint: this often means a missing space after ':' or inconsistent " + "indentation on this line." + ) + elif "could not find expected ':'" in problem_text: + hint = "Hint: a key is likely missing ':' or is mis-indented." + elif "expected " in problem_text: + hint = "Hint: check list/item indentation and unmatched nesting near this line." + + message.append(hint) + return "\n".join(message) + + # ----------Class Definition----------# @@ -1871,8 +1914,9 @@ def __init__(self, test_mode: bool = False) -> None: with open(self.conf_path, "r") as f: try: self.cfg = yaml.load(f, yaml.FullLoader) - except: - logging.critical(__name__ + " : Failed to open config file") + except yaml.YAMLError as error: + formatted = format_yaml_error_message(self.conf_path, error) + logging.critical(__name__ + " : " + formatted) print("Error - See Log") exit() @@ -1882,8 +1926,9 @@ def __init__(self, test_mode: bool = False) -> None: with open(self.sys_path, "r") as f: try: self.sys = yaml.load(f, yaml.FullLoader) - except: - logging.critical(__name__ + " : Failed to open system file") + except yaml.YAMLError as error: + formatted = format_yaml_error_message(self.sys_path, error) + logging.critical(__name__ + " : " + formatted) print("Error - See Log") exit() @@ -1902,10 +1947,22 @@ def yaml_reload(self, test_mode=False) -> Tuple[dict, dict]: if test_mode == False: with open(self.conf_path, "r") as f: - self.cfg = yaml.load(f, yaml.FullLoader) + try: + self.cfg = yaml.load(f, yaml.FullLoader) + except yaml.YAMLError as error: + formatted = format_yaml_error_message(self.conf_path, error) + logging.critical(__name__ + " : " + formatted) + print("Error - See Log") + exit() with open(self.sys_path, "r") as f: - self.sys = yaml.load(f, yaml.FullLoader) + try: + self.sys = yaml.load(f, yaml.FullLoader) + except yaml.YAMLError as error: + formatted = format_yaml_error_message(self.sys_path, error) + logging.critical(__name__ + " : " + formatted) + print("Error - See Log") + exit() return self.cfg, self.sys @@ -1936,8 +1993,9 @@ def __init__(self) -> None: with open(self.parameter_conf, "r") as f: try: self.param = yaml.load(f, yaml.FullLoader) - except: - logging.critical(__name__ + " : Failed to open parameter dictionary") + except yaml.YAMLError as error: + formatted = format_yaml_error_message(self.parameter_conf, error) + logging.critical(__name__ + " : " + formatted) print("Error - See Log") exit() diff --git a/cx_asap/tools/modules/molecule_reconstruction.py b/cx_asap/tools/modules/molecule_reconstruction.py index d8f46f0..5559f2d 100755 --- a/cx_asap/tools/modules/molecule_reconstruction.py +++ b/cx_asap/tools/modules/molecule_reconstruction.py @@ -11,6 +11,10 @@ # ----------Required Modules----------# from system_files.utils import Nice_YAML_Dumper, Config +from system_files.crystal_math import ( + fractional_to_cartesian, + reciprocal_orthonorm_matrix, +) import logging import numpy as np import pathlib @@ -311,22 +315,14 @@ def find_internal_vectors(self) -> None: neutral_coordinates = {} - # The neutral coordinates need to be converted into real space first - - # This is done by multiplying the fractional coordinates by the cell parameters - + # Convert fractional coordinates into Cartesian coordinates using + # full unit-cell orthonormalisation (including non-orthogonal angles). for item in self.neutral_fractional_coordinates: neutral_coordinates[item] = np.array( - [ - [ - self.neutral_fractional_coordinates[item][0] - * self.neutral_cell[0], - self.neutral_fractional_coordinates[item][1] - * self.neutral_cell[1], - self.neutral_fractional_coordinates[item][2] - * self.neutral_cell[2], - ] - ] + fractional_to_cartesian( + self.neutral_fractional_coordinates[item], self.neutral_cell + ), + dtype=float, ) # Vector subtraction to define the bonds based on the construction order @@ -408,18 +404,15 @@ def calculate_fractional_coordinates(self) -> None: self.new_fractional_coordinates = {} self.new_fractional_coordinates[self.starting_atom] = self.starting_coordinates - # First, the internal vectors are converted into fractional coordinates in the new cell + # First, convert internal Cartesian vectors into fractional coordinates + # in the new cell using the reciprocal orthonormal transform. + + M_star = reciprocal_orthonorm_matrix(self.new_cell) for item in self.internal_structure: - new_molecule[item] = [] - new_molecule[item].append( - self.internal_structure[item][0][0] / self.new_cell[0] - ) - new_molecule[item].append( - self.internal_structure[item][0][1] / self.new_cell[1] - ) - new_molecule[item].append( - self.internal_structure[item][0][2] / self.new_cell[2] + new_molecule[item] = np.dot( + np.array(self.internal_structure[item], dtype=float), + M_star.T, ) # This series of loops reconstructs the molecule based on the order previously determined using vector addition @@ -436,9 +429,12 @@ def calculate_fractional_coordinates(self) -> None: for k in self.atomic_order[str(int(i) + 1)]: if k in item.split("-"): self.new_fractional_coordinates[k] = np.add( - self.new_fractional_coordinates[j], + np.array( + self.new_fractional_coordinates[j], + dtype=float, + ), new_molecule[item], - ) + ).tolist() def write_res(self, item: int) -> None: From 4392b8fde53d4fc9ec1bb99914b0a048c480f82f Mon Sep 17 00:00:00 2001 From: Jack Clegg Date: Tue, 14 Jul 2026 11:30:49 +1000 Subject: [PATCH 2/4] fixing an input related issue --- .../modules/centroids.py | 29 +++++++++++++++++-- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/cx_asap/post_refinement_analysis/modules/centroids.py b/cx_asap/post_refinement_analysis/modules/centroids.py index 665e53b..6a8e159 100644 --- a/cx_asap/post_refinement_analysis/modules/centroids.py +++ b/cx_asap/post_refinement_analysis/modules/centroids.py @@ -109,6 +109,7 @@ def _extract_positions( """Gets transformed fractional coordinates for each requested atom.""" atom_names = self._normalise_atom_names(atom_names) + symmetry = self._normalise_symmetry(symmetry) R, t = parse_symm_op(symmetry) if symmetry else (None, None) positions = [] @@ -126,6 +127,24 @@ def _extract_positions( return positions + @staticmethod + def _normalise_symmetry(symmetry: "str | list | None") -> "str | None": + """Converts placeholder symmetry values to None and normalises strings.""" + + if isinstance(symmetry, list): + if len(symmetry) == 0: + return None + symmetry = symmetry[0] + + if symmetry is None: + return None + + value = str(symmetry).strip() + if value in ["", "0", ".", "?"]: + return None + + return value + def calculate_centroid( self, coords: dict, @@ -304,9 +323,13 @@ def point_torsion( logging.warning(__name__ + f" : {error}") return 0.0 - has_symmetry = any( - item is not None for item in [symmetry_1, symmetry_2, symmetry_3, symmetry_4] - ) + symmetry_values = [ + self._normalise_symmetry(symmetry_1), + self._normalise_symmetry(symmetry_2), + self._normalise_symmetry(symmetry_3), + self._normalise_symmetry(symmetry_4), + ] + has_symmetry = any(item is not None for item in symmetry_values) has_centroid_point = any( len(item) > 1 for item in [point_1_list, point_2_list, point_3_list, point_4_list] From 2572b987fdd8de7e407cdc75737ff07f1d39fdd9 Mon Sep 17 00:00:00 2001 From: Jack Clegg Date: Wed, 15 Jul 2026 14:23:34 +1000 Subject: [PATCH 3/4] Make minimum_fraction_of_indexed_spots configurable for pipeline-variable-position --- cx_asap/cxasap.py | 6 ++++++ .../overall_pipelines/variable_position_pipeline.py | 9 ++++++++- cx_asap/system_files/parameter.yaml | 1 + cx_asap/system_files/utils.py | 10 ++++++++-- tests/test_yamls.py | 2 ++ 5 files changed, 25 insertions(+), 3 deletions(-) diff --git a/cx_asap/cxasap.py b/cx_asap/cxasap.py index bd16479..b9dcda6 100644 --- a/cx_asap/cxasap.py +++ b/cx_asap/cxasap.py @@ -402,6 +402,8 @@ def yaml_extraction(heading: str) -> dict: yaml_dict[item] = 8 elif item == "tolerance": yaml_dict[item] = 0.002 + elif item == "minimum_fraction_of_indexed_spots": + yaml_dict[item] = 0.2 elif item == "transformation_matrix": yaml_dict[item] = "1 0 0 0 1 0 0 0 1" elif item == "maximum_cycles": @@ -1001,6 +1003,9 @@ def pipeline_vp(dependencies, files, configure, run): click.echo( " - signal_pixel: enter the values for signal_pixel as a list for XDS processing" ) + click.echo( + " - minimum_fraction_of_indexed_spots: enter the minimum fraction of indexed spots threshold for XDS processing" + ) click.echo( " - structural_analysis_bonds: enter True for bond length analysis, otherwise enter False" ) @@ -1067,6 +1072,7 @@ def pipeline_vp(dependencies, files, configure, run): cfg["atoms_for_rotation_analysis"], cfg["instrument_cif_path"], cfg["total_angle"], + cfg["minimum_fraction_of_indexed_spots"], ) full_vp_analysis.flexible_parameter_loops( diff --git a/cx_asap/overall_pipelines/variable_position_pipeline.py b/cx_asap/overall_pipelines/variable_position_pipeline.py index 06d2cf0..8989468 100755 --- a/cx_asap/overall_pipelines/variable_position_pipeline.py +++ b/cx_asap/overall_pipelines/variable_position_pipeline.py @@ -102,6 +102,7 @@ def setup( MPLA_atoms: str, instrument_cif_path: str, total_angle: int, + minimum_fraction_of_indexed_spots: float, ) -> None: """Organises the directory tree for this experiment @@ -128,6 +129,7 @@ def setup( MPLA_atoms (str): Atoms for mean plane analysis instrument_cif_path (str): full path to the instrument CIF total_angle (int): total wedge angle per experiment + minimum_fraction_of_indexed_spots (float): minimum fraction of indexed spots threshold for XDS indexing """ @@ -139,7 +141,12 @@ def setup( background_reference, instrument_parameters_path, instrument_cif_path ) configure.flexible_setup( - max_processors, neggia_library, space_group_number, MPLA_atoms, total_angle + max_processors, + neggia_library, + space_group_number, + MPLA_atoms, + total_angle, + minimum_fraction_of_indexed_spots, ) self.cfg, self.sys = self.config.yaml_reload() diff --git a/cx_asap/system_files/parameter.yaml b/cx_asap/system_files/parameter.yaml index d3e33a2..bcbb452 100644 --- a/cx_asap/system_files/parameter.yaml +++ b/cx_asap/system_files/parameter.yaml @@ -64,6 +64,7 @@ pipeline-variable-position: - spot_maximum_centroid - signal_pixel - sepmin + - minimum_fraction_of_indexed_spots - mapping_step_size - total_angle - calculate_centroid_distance diff --git a/cx_asap/system_files/utils.py b/cx_asap/system_files/utils.py index c33ca82..117ee8a 100755 --- a/cx_asap/system_files/utils.py +++ b/cx_asap/system_files/utils.py @@ -143,6 +143,7 @@ def flexible_setup( space_group_number: int, MPLA_atoms: str, total_angle: int, + minimum_fraction_of_indexed_spots: float = 0.2, ) -> None: """Primarily sets up the XDS.INP file for a flexible crystal experiment @@ -155,6 +156,7 @@ def flexible_setup( MPLA_atoms (str): atoms for MPLA command to be written into reference .ins/.res file total_angle (int): total wedge angle measured in the experiment + minimum_fraction_of_indexed_spots (float): minimum fraction of indexed spots threshold for XDS indexing """ @@ -236,7 +238,7 @@ def flexible_setup( self.XDS.change( self.sys["XDS_inp_organised"], "MINIMUM_FRACTION_OF_INDEXED_SPOTS", - 0.2, + minimum_fraction_of_indexed_spots, ) flag1 += 1 elif "SEPMIN" in line: @@ -271,7 +273,11 @@ def flexible_setup( with open(self.sys["XDS_inp_organised"], "a") as in_file: if flag1 == 0: - in_file.write(" MINIMUM_FRACTION_OF_INDEXED_SPOTS= 0.2\n") + in_file.write( + " MINIMUM_FRACTION_OF_INDEXED_SPOTS= " + + str(minimum_fraction_of_indexed_spots) + + "\n" + ) if flag2 == 0: in_file.write(" SEPMIN= 7\n") if flag3 == 0: diff --git a/tests/test_yamls.py b/tests/test_yamls.py index 1ab2d4c..7bf8c54 100644 --- a/tests/test_yamls.py +++ b/tests/test_yamls.py @@ -69,6 +69,7 @@ def setUp(self): "spot_maximum_centroid", "signal_pixel", "sepmin", + "minimum_fraction_of_indexed_spots", "mapping_step_size", "total_angle", "calculate_centroid_distance", @@ -439,6 +440,7 @@ def setUp(self): "spot_maximum_centroid", "signal_pixel", "sepmin", + "minimum_fraction_of_indexed_spots", "mapping_step_size", "total_angle", "calculate_centroid_distance", From a7deed5c8a99fa52a5fbca47d9e12f111c077ad7 Mon Sep 17 00:00:00 2001 From: Jack Clegg Date: Thu, 30 Jul 2026 10:20:10 +1000 Subject: [PATCH 4/4] updating cif-analysis to be consistent with new points changes --- cx_asap/cxasap.py | 227 +----- .../variable_position_pipeline.py | 14 +- .../variable_temperature_pipeline.py | 14 +- .../modules/cif_geometry.py | 93 +-- .../modules/points.py | 685 ++++++++++++++++++ .../pipelines/points_pipeline.py | 280 +++++++ cx_asap/system_files/parameter.yaml | 66 +- tests/test_yamls.py | 134 +++- 8 files changed, 1161 insertions(+), 352 deletions(-) create mode 100644 cx_asap/post_refinement_analysis/modules/points.py create mode 100644 cx_asap/post_refinement_analysis/pipelines/points_pipeline.py diff --git a/cx_asap/cxasap.py b/cx_asap/cxasap.py index b9dcda6..4e2786b 100644 --- a/cx_asap/cxasap.py +++ b/cx_asap/cxasap.py @@ -133,10 +133,10 @@ from post_refinement_analysis.modules.rotation_planes import Rotation from post_refinement_analysis.modules.structural_analysis import Structural_Analysis from post_refinement_analysis.modules.ADP_analysis import ADP_analysis -from post_refinement_analysis.modules.centroids import Centroids +from post_refinement_analysis.modules.points import PointGeometryEngine from post_refinement_analysis.modules.cif_analysis import CIF_Analysis from post_refinement_analysis.pipelines.rotation_pipeline import Rotation_Pipeline -from post_refinement_analysis.pipelines.centroids_pipeline import Centroids_Pipeline +from post_refinement_analysis.pipelines.points_pipeline import PointsPipeline from post_refinement_analysis.pipelines.cif_analysis_pipeline import CIF_Analysis_Pipeline from post_refinement_analysis.pipelines.variable_cif_parameter import ( Variable_Analysis_Pipeline, @@ -456,8 +456,8 @@ def configuration_check(heading: str) -> Tuple[bool, dict]: "beta_gradient", "c_gradient", "gamma_gradient", - "centroid_1_symmetry", - "centroid_2_symmetry", + "point_group_1_symmetry", + "point_group_2_symmetry", ] if heading == "pipeline-AS-Brute-individual": @@ -1025,19 +1025,19 @@ def pipeline_vp(dependencies, files, configure, run): " - wedge_angles: enter the wedge angles as a list for XDS processing" ) click.echo( - " - calculate_centroid_distance: enter True for centroid analysis between two atom-group centroids from the .lst files, otherwise enter False" + " - calculate_point_group_distance: enter True for point-group analysis between two atom groups from the .lst files, otherwise enter False" ) click.echo( - " - centroid_1_atoms: list of atom labels for the first centroid group (only needed if calculate_centroid_distance is true)" + " - point_group_1_atoms: list of atom labels for the first point group (only needed if calculate_point_group_distance is true)" ) click.echo( - " - centroid_2_atoms: list of atom labels for the second centroid group (only needed if calculate_centroid_distance is true)" + " - point_group_2_atoms: list of atom labels for the second point group (only needed if calculate_point_group_distance is true)" ) click.echo( - " - centroid_1_symmetry: optional symmetry operation for centroid 1 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" + " - point_group_1_symmetry: optional symmetry operation for point group 1 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" ) click.echo( - " - centroid_2_symmetry: optional symmetry operation for centroid 2 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" + " - point_group_2_symmetry: optional symmetry operation for point group 2 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" ) fields = yaml_extraction("pipeline-variable-position") @@ -1122,15 +1122,15 @@ def pipeline_vp(dependencies, files, configure, run): cfg["wedge_angles"], cfg["reference_plane"], ) - if cfg.get("calculate_centroid_distance", False): - centroid_analysis = Centroids_Pipeline() - centroid_analysis.centroid_distance_analysis( + if cfg.get("calculate_point_group_distance", False): + centroid_analysis = PointsPipeline() + centroid_analysis.point_group_distance_analysis( full_vp_analysis.sys["current_results_path"], - cfg["centroid_1_atoms"], - cfg["centroid_2_atoms"], + cfg["point_group_1_atoms"], + cfg["point_group_2_atoms"], full_vp_analysis.sys["current_results_path"], - symmetry_1=cfg.get("centroid_1_symmetry") or None, - symmetry_2=cfg.get("centroid_2_symmetry") or None, + symmetry_1=cfg.get("point_group_1_symmetry") or None, + symmetry_2=cfg.get("point_group_2_symmetry") or None, ) copy_logs(full_vp_analysis.sys["current_results_path"]) @@ -3491,173 +3491,6 @@ def pipeline_rotation_planes(dependencies, files, configure, run): click.echo("Please select an option. To view options, add --help") -######----- Module Centroids ------##### - - -@click.command( - "module-centroids", - short_help="calculate centroid distance between two atom groups", -) -@click.option("--dependencies", is_flag=True, help="view the software dependencies") -@click.option("--files", is_flag=True, help="view the required input files") -@click.option("--configure", is_flag=True, help="generate your conf.yaml file") -@click.option("--run", is_flag=True, help="run the code!") -def module_centroids(dependencies, files, configure, run): - """For a single dataset, calculate the distance between two atom-group centroids - and/or the SHELXL inter-plane angle from a .lst file. - """ - if dependencies: - click.echo("\nYou do not require any additional software in your path!\n") - elif files: - click.echo("\nYou require the below files:") - click.echo(" - a .lst file output after refinement in SHELXL") - click.echo(" - for inter-plane angle, the .lst must contain two MPLA commands") - click.echo("\nThis file can be located anywhere ") - elif configure: - click.echo("\nWriting a file called conf.yaml in the cx_asap folder...\n") - click.echo("You will need to fill out the parameters.") - click.echo("Descriptions are listed below:") - click.echo( - " - lst_file_location: enter the full path to your lst file for analysis" - ) - click.echo( - " - centroid_1_atoms: list of atom labels for the first centroid group" - ) - click.echo( - " - centroid_2_atoms: list of atom labels for the second centroid group" - ) - click.echo( - " - centroid_1_symmetry: optional symmetry operation for centroid 1 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" - ) - click.echo( - " - centroid_2_symmetry: optional symmetry operation for centroid 2 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" - ) - - fields = yaml_extraction("module-centroids") - yaml_creation(fields) - - elif run: - click.echo("\nChecking to see if experiment configured....\n") - - check, cfg = configuration_check("module-centroids") - - if check == False: - click.echo("Make sure you fill in the configuration file!") - click.echo( - "If you last ran a different code, make sure you reconfigure for the new script!" - ) - click.echo("Re-run configuration for description of each parameter\n") - else: - click.echo("READY TO RUN SCRIPT!\n") - reset_logs() - results_dir = pathlib.Path(cfg["lst_file_location"]).parent - centroid_analysis = Centroids() - centroid_analysis.analyse_centroid_distance( - cfg["lst_file_location"], - 1, - results_dir, - cfg["centroid_1_atoms"], - cfg["centroid_2_atoms"], - symmetry_1=cfg.get("centroid_1_symmetry") or None, - symmetry_2=cfg.get("centroid_2_symmetry") or None, - ) - - copy_logs(results_dir) - - output_message() - - else: - click.echo("Please select an option. To view options, add --help") - - -#####------ Pipeline Centroids -----####### - - -@click.command( - "pipeline-centroids", - short_help="calculate centroid distances for multiple datasets", -) -@click.option("--dependencies", is_flag=True, help="view the software dependencies") -@click.option("--files", is_flag=True, help="view the required input files") -@click.option("--configure", is_flag=True, help="generate your conf.yaml file") -@click.option("--run", is_flag=True, help="run the code!") -def pipeline_centroids(dependencies, files, configure, run): - """For a series of datasets, calculate centroid distances and/or SHELXL - inter-plane angles from .lst files across multiple folders. - """ - if dependencies: - click.echo("\nYou do not require any additional software in your path!\n") - elif files: - click.echo("\nYou require the below files:") - click.echo( - " - a series of .lst files in separate folders contained in a single parent folder" - ) - click.echo(" - for inter-plane angle, each .lst must contain two MPLA commands") - click.echo("\nThis parent folder can be located anywhere ") - click.echo( - "Results will be written to a numbered folder inside Centroid_Analysis within this parent folder" - ) - elif configure: - click.echo("\nWriting a file called conf.yaml in the cx_asap folder...\n") - click.echo("You will need to fill out the parameters.") - click.echo("Descriptions are listed below:") - click.echo( - " - experiment_location: full path to the parent folder containing a series of folders with .lst files inside" - ) - click.echo( - " - output: results will be written to a numbered folder inside Centroid_Analysis in the experiment_location folder" - ) - click.echo( - " - centroid_1_atoms: list of atom labels for the first centroid group" - ) - click.echo( - " - centroid_2_atoms: list of atom labels for the second centroid group" - ) - click.echo( - " - centroid_1_symmetry: optional symmetry operation for centroid 1 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" - ) - click.echo( - " - centroid_2_symmetry: optional symmetry operation for centroid 2 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" - ) - - fields = yaml_extraction("pipeline-centroids") - yaml_creation(fields) - - elif run: - click.echo("\nChecking to see if experiment configured....\n") - - check, cfg = configuration_check("pipeline-centroids") - - if check == False: - click.echo("Make sure you fill in the configuration file!") - click.echo( - "If you last ran a different code, make sure you reconfigure for the new script!" - ) - click.echo("Re-run configuration for description of each parameter\n") - else: - click.echo("READY TO RUN SCRIPT!\n") - reset_logs() - multi_centroid = Centroids_Pipeline() - results_dir = multi_centroid.create_numbered_results_directory( - cfg["experiment_location"], "Centroid_Analysis" - ) - multi_centroid.centroid_distance_analysis( - cfg["experiment_location"], - cfg["centroid_1_atoms"], - cfg["centroid_2_atoms"], - results_dir, - symmetry_1=cfg.get("centroid_1_symmetry") or None, - symmetry_2=cfg.get("centroid_2_symmetry") or None, - ) - - copy_logs(results_dir) - - output_message() - - else: - click.echo("Please select an option. To view options, add --help") - - #####------ Module Point Geometry -----####### @@ -3778,7 +3611,7 @@ def module_point_geometry(dependencies, files, configure, run): return results_dir = pathlib.Path(cfg["lst_file_location"]).parent - centroid_analysis = Centroids() + centroid_analysis = PointGeometryEngine() centroid_analysis.analyse_point_geometry( cfg["lst_file_location"], 1, @@ -3925,7 +3758,7 @@ def pipeline_point_geometry(dependencies, files, configure, run): output_message() return - multi_geometry = Centroids_Pipeline() + multi_geometry = PointsPipeline() results_dir = multi_geometry.create_numbered_results_directory( cfg["experiment_location"], "Geometry_Analysis" ) @@ -4565,19 +4398,19 @@ def pipeline_position_analysis(dependencies, files, configure, run): " - wedge_angles: enter the wedge angles as a list for XDS processing" ) click.echo( - " - calculate_centroid_distance: enter True for centroid analysis between two atom-group centroids from the .lst files, otherwise enter False" + " - calculate_point_group_distance: enter True for point-group analysis between two atom groups from the .lst files, otherwise enter False" ) click.echo( - " - centroid_1_atoms: list of atom labels for the first centroid group (only needed if calculate_centroid_distance is true)" + " - point_group_1_atoms: list of atom labels for the first point group (only needed if calculate_point_group_distance is true)" ) click.echo( - " - centroid_2_atoms: list of atom labels for the second centroid group (only needed if calculate_centroid_distance is true)" + " - point_group_2_atoms: list of atom labels for the second point group (only needed if calculate_point_group_distance is true)" ) click.echo( - " - centroid_1_symmetry: optional symmetry operation for centroid 1 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" + " - point_group_1_symmetry: optional symmetry operation for point group 1 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" ) click.echo( - " - centroid_2_symmetry: optional symmetry operation for centroid 2 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" + " - point_group_2_symmetry: optional symmetry operation for point group 2 atoms, e.g. -x+1/2, y+1/2, -z+1/2 (leave blank if not required)" ) fields = yaml_extraction("pipeline-position-analysis") @@ -4617,15 +4450,15 @@ def pipeline_position_analysis(dependencies, files, configure, run): cfg["structural_analysis_hbonds"], cfg["ADP_analysis"], ) - if cfg.get("calculate_centroid_distance", False): - centroid_analysis = Centroids_Pipeline() - centroid_analysis.centroid_distance_analysis( + if cfg.get("calculate_point_group_distance", False): + centroid_analysis = PointsPipeline() + centroid_analysis.point_group_distance_analysis( cfg["experiment_location"], - cfg["centroid_1_atoms"], - cfg["centroid_2_atoms"], + cfg["point_group_1_atoms"], + cfg["point_group_2_atoms"], cfg["experiment_location"], - symmetry_1=cfg.get("centroid_1_symmetry") or None, - symmetry_2=cfg.get("centroid_2_symmetry") or None, + symmetry_1=cfg.get("point_group_1_symmetry") or None, + symmetry_2=cfg.get("point_group_2_symmetry") or None, ) copy_logs(cfg["experiment_location"]) diff --git a/cx_asap/overall_pipelines/variable_position_pipeline.py b/cx_asap/overall_pipelines/variable_position_pipeline.py index 8989468..fe2508c 100755 --- a/cx_asap/overall_pipelines/variable_position_pipeline.py +++ b/cx_asap/overall_pipelines/variable_position_pipeline.py @@ -23,7 +23,7 @@ from data_refinement.pipelines.refine_pipeline import Refinement_Pipeline from cif_validation.pipelines.cif_pipeline import CIF_Compile_Pipeline from post_refinement_analysis.pipelines.rotation_pipeline import Rotation_Pipeline -from post_refinement_analysis.pipelines.centroids_pipeline import Centroids_Pipeline +from post_refinement_analysis.pipelines.points_pipeline import PointsPipeline from post_refinement_analysis.pipelines.variable_position_analysis import ( VP_Analysis_Pipeline, ) @@ -530,14 +530,14 @@ def process( ) rotation = Rotation_Pipeline() rotation.analysis(location, reference_plane, graph_output_location) - centroids = Centroids_Pipeline() - centroids.centroid_distance_analysis( + centroids = PointsPipeline() + centroids.point_group_distance_analysis( location, - self.cfg["centroid_1_atoms"], - self.cfg["centroid_2_atoms"], + self.cfg["point_group_1_atoms"], + self.cfg["point_group_2_atoms"], graph_output_location, - symmetry_1=self.cfg.get("centroid_1_symmetry"), - symmetry_2=self.cfg.get("centroid_2_symmetry"), + symmetry_1=self.cfg.get("point_group_1_symmetry"), + symmetry_2=self.cfg.get("point_group_2_symmetry"), ) cif = CIF_Compile_Pipeline() cif.configure( diff --git a/cx_asap/overall_pipelines/variable_temperature_pipeline.py b/cx_asap/overall_pipelines/variable_temperature_pipeline.py index cc58ac0..1bd3adf 100755 --- a/cx_asap/overall_pipelines/variable_temperature_pipeline.py +++ b/cx_asap/overall_pipelines/variable_temperature_pipeline.py @@ -18,7 +18,7 @@ VT_Analysis_Pipeline, ) from post_refinement_analysis.pipelines.rotation_pipeline import Rotation_Pipeline -from post_refinement_analysis.pipelines.centroids_pipeline import Centroids_Pipeline +from post_refinement_analysis.pipelines.points_pipeline import PointsPipeline import yaml import os import logging @@ -142,14 +142,14 @@ def process( ) rotation = Rotation_Pipeline() rotation.analysis(location, reference_plane, graph_output_location) - centroids = Centroids_Pipeline() - centroids.centroid_distance_analysis( + centroids = PointsPipeline() + centroids.point_group_distance_analysis( location, - self.cfg["centroid_1_atoms"], - self.cfg["centroid_2_atoms"], + self.cfg["point_group_1_atoms"], + self.cfg["point_group_2_atoms"], graph_output_location, - symmetry_1=self.cfg.get("centroid_1_symmetry"), - symmetry_2=self.cfg.get("centroid_2_symmetry"), + symmetry_1=self.cfg.get("point_group_1_symmetry"), + symmetry_2=self.cfg.get("point_group_2_symmetry"), ) def analyse( diff --git a/cx_asap/post_refinement_analysis/modules/cif_geometry.py b/cx_asap/post_refinement_analysis/modules/cif_geometry.py index 8ba4865..aab081b 100644 --- a/cx_asap/post_refinement_analysis/modules/cif_geometry.py +++ b/cx_asap/post_refinement_analysis/modules/cif_geometry.py @@ -9,7 +9,8 @@ ################################################################################################### from CifFile import ReadCif -from post_refinement_analysis.modules.centroids import Centroids +from post_refinement_analysis.modules.cif_read import CIF_Read +from post_refinement_analysis.modules.points import PointGeometryEngine from system_files.crystal_math import ( fractional_to_cartesian, best_fit_plane_normal, @@ -26,10 +27,10 @@ class CIF_Geometry: """Performs symmetry-aware point and plane analysis directly from CIF files.""" def __init__(self, test_mode: bool = False) -> None: - """Initialise CIF geometry engine and centroid helper.""" + """Initialise CIF geometry engine and point-geometry helper.""" self.test_mode = test_mode - self.point_engine = Centroids(self.test_mode) + self.point_engine = PointGeometryEngine(self.test_mode) @staticmethod def _parse_float(raw) -> float: @@ -47,72 +48,6 @@ def _parse_float(raw) -> float: return float(raw_str) - @staticmethod - def _read_cif_files(cif_location: str) -> list: - """Discover CIF inputs from a file path or root-first folder layout.""" - - base = pathlib.Path(cif_location) - - if base.is_file(): - if base.suffix.lower() == ".cif": - return [base.resolve()] - return [] - - if not base.exists() or not base.is_dir(): - return [] - - def _is_results_folder(folder_name: str) -> bool: - name = folder_name.strip().lower() - blocked_exact = { - "cif_analysis", - "geometry_analysis", - "refinement_statistics", - "results", - "analysis", - "ref", - "failed_autoprocessing", - } - if name in blocked_exact: - return True - if name.startswith("_"): - return True - return False - - root_files = [item for item in sorted(base.glob("*.cif")) if item.is_file()] - - files = [] - - # One folder level down only - for child in sorted(base.iterdir()): - if not child.is_dir(): - continue - if _is_results_folder(child.name): - continue - files.extend( - [item for item in sorted(child.glob("*.cif")) if item.is_file()] - ) - - if len(root_files) > 0: - if len(files) > 0: - logging.warning( - __name__ - + " : Mixed CIF layout detected (root-level and nested CIFs). " - + "Using root-level CIFs only; nested CIFs will be ignored." - ) - return [item.resolve() for item in root_files] - - # Deduplicate while preserving deterministic order - seen = set() - unique_files = [] - for item in files: - key = str(item.resolve()) - if key in seen: - continue - seen.add(key) - unique_files.append(item) - - return unique_files - def _extract_structure_data(self, block) -> dict: """Extract unit-cell and fractional atom coordinates from one CIF block.""" @@ -166,8 +101,10 @@ def _append_rows(csv_name: str, rows: list, results_directory: str) -> None: df.to_csv(output_path, index=None) @staticmethod - def _definition_has_centroid(definition: dict, keys: list, engine: Centroids) -> bool: - """Check whether any definition key uses centroid-style atom syntax.""" + def _definition_uses_mercury_rounding( + definition: dict, keys: list, engine: PointGeometryEngine + ) -> bool: + """Check whether any definition key should also be written in Mercury-rounded form.""" for key in keys: if engine._point_uses_centroid(definition.get(key, [])): @@ -243,8 +180,8 @@ def run( For each CIF block, this writes CSV outputs for any configured distance, angle, torsion, point-plane, rotation-plane, and optional interplane - calculations. Mercury-style rounded-centroid outputs are generated when - ``mercury_output`` is enabled and centroid definitions are present. + calculations. Mercury-style rounded outputs are generated when + ``mercury_output`` is enabled and definitions require companion rounded values. """ point_geometry_distances = point_geometry_distances or [] @@ -254,7 +191,7 @@ def run( rotation_plane_definitions = rotation_plane_definitions or [] mean_plane_definitions = mean_plane_definitions or [] - files = self._read_cif_files(cif_location) + files = CIF_Read._read_cif_files(cif_location) if len(files) == 0: logging.warning(__name__ + " : No .cif files found for CIF geometry analysis") return @@ -325,7 +262,7 @@ def run( definition.get("point_1_symmetry") or None, definition.get("point_2_symmetry") or None, ) - if mercury_output and self._definition_has_centroid( + if mercury_output and self._definition_uses_mercury_rounding( definition, ["point_1_atoms", "point_2_atoms"], self.point_engine, @@ -352,7 +289,7 @@ def run( definition.get("point_2_symmetry") or None, definition.get("point_3_symmetry") or None, ) - if mercury_output and self._definition_has_centroid( + if mercury_output and self._definition_uses_mercury_rounding( definition, ["point_1_atoms", "point_2_atoms", "point_3_atoms"], self.point_engine, @@ -383,7 +320,7 @@ def run( definition.get("point_3_symmetry") or None, definition.get("point_4_symmetry") or None, ) - if mercury_output and self._definition_has_centroid( + if mercury_output and self._definition_uses_mercury_rounding( definition, [ "point_1_atoms", @@ -417,7 +354,7 @@ def run( definition.get("point_symmetry") or None, definition.get("plane_symmetry") or None, ) - if mercury_output and self._definition_has_centroid( + if mercury_output and self._definition_uses_mercury_rounding( definition, ["point_atoms"], self.point_engine, diff --git a/cx_asap/post_refinement_analysis/modules/points.py b/cx_asap/post_refinement_analysis/modules/points.py new file mode 100644 index 0000000..e9e8c93 --- /dev/null +++ b/cx_asap/post_refinement_analysis/modules/points.py @@ -0,0 +1,685 @@ +#!/usr/bin/env python3 + +################################################################################################### +# -----------------------------------------CX-ASAP: points----------------------------------------# +# ---Authors: Amy J. Thompson, Kate M. Smith, Daniel J. Eriksson, Jack K. Clegg & Jason R. Price---# +# -----------------------------------Python Implementation by AJT----------------------------------# +# -----------------------------------Project Design by JRP and JKC---------------------------------# +# --------------------------------Valuable Coding Support by KMS & DJE-----------------------------# +################################################################################################### + +# ----------Required Modules----------# + +from system_files.utils import Config +from system_files.crystal_math import ( + parse_symm_op, + fractional_to_cartesian, + distance_between_points, + angle_between_points, + torsion_between_points, + point_to_plane_distance, +) +from post_refinement_analysis.modules.lst_read import LST_Read +import pathlib +import os +import logging +import numpy as np +import pandas as pd + +# ----------Class Definition----------# + + +class PointGeometryEngine: + """Calculates centroid positions and inter-centroid distances from .lst files. + + Atom fractional coordinates are parsed from the embedded .res block of the + SHELXL .lst file. Distances are computed in Cartesian coordinates using the + unit cell orthonormalisation matrix. + + This class can also extract the inter-plane angle reported by SHELXL when + two MPLA commands are present. + """ + + def __init__(self, test_mode: bool = False) -> None: + """Initialises the class. + + Args: + test_mode (bool): if True, skips conf.yaml loading + """ + + self.test_mode = test_mode + + config = Config(self.test_mode) + self.cfg = config.cfg + self.sys = config.sys + self.conf_path = config.conf_path + self.sys_path = config.sys_path + + self.lst_reader = LST_Read(self.test_mode) + + def grab_cell(self, file_name: str) -> None: + """Reads the unit cell parameters from the CELL line of a .lst file. + + Args: + file_name (str): full path to the .lst file + """ + + self.bad_flag = False + + with open(file_name, "rt") as f: + split_line = [] + for line in f: + if line.startswith(" CELL") or line.startswith("CELL"): + split_line = line.split() + break + + if len(split_line) >= 8: + self.cell_params = [float(split_line[i]) for i in range(2, 8)] + else: + logging.warning( + __name__ + " : Could not read CELL line from " + str(file_name) + ) + self.bad_flag = True + + def _normalise_atom_names( + self, atom_names: "str | list[str] | tuple | set | np.ndarray" + ) -> list: + """Normalises atom label input into a list of strings.""" + + if isinstance(atom_names, str): + atom_names = atom_names.split() + elif isinstance(atom_names, (list, tuple, set, np.ndarray)): + atom_names = list(atom_names) + else: + message = ( + "atom_names must be a list/tuple/set/ndarray of labels or " + "a whitespace-delimited string" + ) + logging.error(__name__ + f" : {message}. Got {type(atom_names).__name__}") + raise TypeError(message) + + return [str(item) for item in atom_names if str(item).strip()] + + def _extract_positions( + self, + coords: dict, + atom_names: "str | list[str] | tuple | set | np.ndarray", + symmetry: str = None, + ) -> list: + """Gets transformed fractional coordinates for each requested atom.""" + + atom_names = self._normalise_atom_names(atom_names) + R, t = parse_symm_op(symmetry) if symmetry else (None, None) + + positions = [] + for name in atom_names: + key = name.upper() + if key in coords: + pos = np.array(coords[key], dtype=float) + if R is not None: + pos = np.dot(R.T, pos) + t + positions.append(pos) + else: + logging.warning( + __name__ + f" : Atom {name} not found in coordinate list" + ) + + return positions + + def calculate_point_group_center( + self, + coords: dict, + atom_names: "str | list[str] | tuple | set | np.ndarray", + symmetry: str = None, + ) -> "np.ndarray | None": + """Calculates the centroid (mean fractional position) of a group of atoms. + + If a symmetry operation string is provided, each atom's fractional + coordinates are transformed by that operation before averaging. + + Args: + coords (dict): fractional coordinates from LST_Read.extract_atom_coordinates() + atom_names (list): atom labels to include in the centroid calculation + symmetry (str): optional SHELXL/CIF symmetry operation string, + e.g. "-x+1/2, y+1/2, -z+1/2" + + Returns: + centroid (np.ndarray): 1D array [x, y, z] in fractional coordinates, + or None if no atoms were found + """ + + positions = self._extract_positions(coords, atom_names, symmetry) + + if not positions: + logging.critical( + __name__ + " : No valid atoms found for centroid calculation" + ) + return None + + return np.mean(positions, axis=0) + + def calculate_point( + self, + coords: dict, + atom_names: "str | list[str] | tuple | set | np.ndarray", + symmetry: str = None, + group_center_round_dp: int = None, + ) -> "np.ndarray | None": + """Calculates a single point from atom input. + + One atom label resolves to that atom position. + Multiple atom labels resolve to the centroid of those atoms. + """ + + atom_names = self._normalise_atom_names(atom_names) + if not atom_names: + logging.critical(__name__ + " : No atoms supplied for point calculation") + return None + + if len(atom_names) == 1: + positions = self._extract_positions(coords, atom_names, symmetry) + if not positions: + logging.critical( + __name__ + " : Could not resolve atom for point calculation" + ) + return None + return positions[0] + + centroid = self.calculate_point_group_center(coords, atom_names, symmetry) + if centroid is None: + return None + + if group_center_round_dp is not None: + return np.round(centroid, group_center_round_dp) + + return centroid + + def _point_uses_group_center( + self, atom_names: "str | list[str] | tuple | set | np.ndarray" + ) -> bool: + """Returns True when a point definition resolves via centroid averaging.""" + + return len(self._normalise_atom_names(atom_names)) > 1 + + def _fractional_to_cartesian(self, point_frac: "np.ndarray | list") -> np.ndarray: + """Converts one fractional point into Cartesian coordinates.""" + + return fractional_to_cartesian(point_frac, self.cell_params) + + def point_distance( + self, + coords: dict, + point_1_atoms: "str | list[str] | tuple | set | np.ndarray", + point_2_atoms: "str | list[str] | tuple | set | np.ndarray", + symmetry_1: str = None, + symmetry_2: str = None, + group_center_round_dp: int = None, + ) -> float: + """Calculates Cartesian distance between two points.""" + + p1 = self.calculate_point(coords, point_1_atoms, symmetry_1, group_center_round_dp) + p2 = self.calculate_point(coords, point_2_atoms, symmetry_2, group_center_round_dp) + + if p1 is None or p2 is None: + return 0.0 + + cart1 = self._fractional_to_cartesian(p1) + cart2 = self._fractional_to_cartesian(p2) + return distance_between_points(cart1, cart2) + + def point_angle( + self, + coords: dict, + point_1_atoms: "str | list[str] | tuple | set | np.ndarray", + point_2_atoms: "str | list[str] | tuple | set | np.ndarray", + point_3_atoms: "str | list[str] | tuple | set | np.ndarray", + symmetry_1: str = None, + symmetry_2: str = None, + symmetry_3: str = None, + group_center_round_dp: int = None, + ) -> float: + """Calculates angle (degrees) formed by points 1-2-3 at point 2.""" + + p1 = self.calculate_point(coords, point_1_atoms, symmetry_1, group_center_round_dp) + p2 = self.calculate_point(coords, point_2_atoms, symmetry_2, group_center_round_dp) + p3 = self.calculate_point(coords, point_3_atoms, symmetry_3, group_center_round_dp) + + if p1 is None or p2 is None or p3 is None: + return 0.0 + + c1 = self._fractional_to_cartesian(p1) + c2 = self._fractional_to_cartesian(p2) + c3 = self._fractional_to_cartesian(p3) + + try: + return angle_between_points(c1, c2, c3) + except ValueError as error: + logging.warning(__name__ + f" : {error}") + return 0.0 + + def point_torsion( + self, + coords: dict, + point_1_atoms: "str | list[str] | tuple | set | np.ndarray", + point_2_atoms: "str | list[str] | tuple | set | np.ndarray", + point_3_atoms: "str | list[str] | tuple | set | np.ndarray", + point_4_atoms: "str | list[str] | tuple | set | np.ndarray", + symmetry_1: str = None, + symmetry_2: str = None, + symmetry_3: str = None, + symmetry_4: str = None, + group_center_round_dp: int = None, + ) -> float: + """Calculates torsion angle (degrees) for points 1-2-3-4.""" + + point_1_list = self._normalise_atom_names(point_1_atoms) + point_2_list = self._normalise_atom_names(point_2_atoms) + point_3_list = self._normalise_atom_names(point_3_atoms) + point_4_list = self._normalise_atom_names(point_4_atoms) + + p1 = self.calculate_point( + coords, point_1_list, symmetry_1, group_center_round_dp + ) + p2 = self.calculate_point( + coords, point_2_list, symmetry_2, group_center_round_dp + ) + p3 = self.calculate_point( + coords, point_3_list, symmetry_3, group_center_round_dp + ) + p4 = self.calculate_point( + coords, point_4_list, symmetry_4, group_center_round_dp + ) + + if p1 is None or p2 is None or p3 is None or p4 is None: + return 0.0 + + c1 = self._fractional_to_cartesian(p1) + c2 = self._fractional_to_cartesian(p2) + c3 = self._fractional_to_cartesian(p3) + c4 = self._fractional_to_cartesian(p4) + + try: + torsion_value = torsion_between_points(c1, c2, c3, c4) + except ValueError as error: + logging.warning(__name__ + f" : {error}") + return 0.0 + + has_symmetry = any( + item is not None for item in [symmetry_1, symmetry_2, symmetry_3, symmetry_4] + ) + has_centroid_point = any( + len(item) > 1 + for item in [point_1_list, point_2_list, point_3_list, point_4_list] + ) + + if has_symmetry and has_centroid_point and abs(torsion_value) < 90.0: + return torsion_value - 180.0 if torsion_value > 0.0 else torsion_value + 180.0 + + return torsion_value + + def point_plane_distance( + self, + coords: dict, + point_atoms: "str | list[str] | tuple | set | np.ndarray", + plane_atoms: "str | list[str] | tuple | set | np.ndarray", + point_symmetry: str = None, + plane_symmetry: str = None, + group_center_round_dp: int = None, + ) -> float: + """Calculates absolute distance from a point to a best-fit plane.""" + + point_frac = self.calculate_point( + coords, point_atoms, point_symmetry, group_center_round_dp + ) + if point_frac is None: + return 0.0 + + plane_positions = self._extract_positions(coords, plane_atoms, plane_symmetry) + if len(plane_positions) < 3: + logging.warning(__name__ + " : Need at least 3 atoms to define a plane") + return 0.0 + + point_cart = self._fractional_to_cartesian(point_frac) + plane_cart = np.array( + [self._fractional_to_cartesian(item) for item in plane_positions], + dtype=float, + ) + + try: + return point_to_plane_distance(point_cart, plane_cart) + except ValueError as error: + logging.warning(__name__ + f" : {error}") + return 0.0 + + def _append_measurements_csv( + self, + csv_name: str, + structure_number: int, + values: "dict[str, float]", + results_path: str, + ) -> None: + """Appends one structure worth of labelled measurements to a CSV.""" + + if not values: + return + + row = {"Structure": [structure_number]} + for label, value in values.items(): + row[label] = [value] + + df = pd.DataFrame(row) + + os.chdir(results_path) + try: + old_data = pd.read_csv(csv_name) + except FileNotFoundError: + df.to_csv(csv_name, index=None) + else: + new_df = pd.concat([old_data, df]) + new_df.to_csv(csv_name, index=None) + + def analyse_point_geometry( + self, + lst_name: str, + structure_number: int, + results_path: str, + distance_definitions: "list[dict]" = None, + angle_definitions: "list[dict]" = None, + torsion_definitions: "list[dict]" = None, + plane_distance_definitions: "list[dict]" = None, + mercury_output: bool = False, + ) -> None: + """Calculates configured point-geometry values and appends CSV outputs.""" + + if lst_name == "": + logging.info(__name__ + " : Refinement failed, no structure to analyse") + return + + self.grab_cell(pathlib.Path(lst_name)) + + if self.bad_flag: + return + + data = self.lst_reader.read(lst_name) + coords = self.lst_reader.extract_atom_coordinates(data) + + distance_values = {} + mercury_distance_values = {} + angle_values = {} + mercury_angle_values = {} + torsion_values = {} + mercury_torsion_values = {} + plane_distance_values = {} + mercury_plane_distance_values = {} + + for index, definition in enumerate(distance_definitions or []): + if not isinstance(definition, dict): + continue + label = definition.get("label", f"Distance_{index + 1}") + distance_values[label] = self.point_distance( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + ) + if mercury_output and any( + self._point_uses_group_center(definition.get(key, [])) + for key in ["point_1_atoms", "point_2_atoms"] + ): + mercury_distance_values[label] = self.point_distance( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + group_center_round_dp=3, + ) + + for index, definition in enumerate(angle_definitions or []): + if not isinstance(definition, dict): + continue + label = definition.get("label", f"Angle_{index + 1}") + angle_values[label] = self.point_angle( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_3_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + definition.get("point_3_symmetry") or None, + ) + if mercury_output and any( + self._point_uses_group_center(definition.get(key, [])) + for key in ["point_1_atoms", "point_2_atoms", "point_3_atoms"] + ): + mercury_angle_values[label] = self.point_angle( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_3_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + definition.get("point_3_symmetry") or None, + group_center_round_dp=3, + ) + + for index, definition in enumerate(torsion_definitions or []): + if not isinstance(definition, dict): + continue + label = definition.get("label", f"Torsion_{index + 1}") + torsion_values[label] = self.point_torsion( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_3_atoms", []), + definition.get("point_4_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + definition.get("point_3_symmetry") or None, + definition.get("point_4_symmetry") or None, + ) + if mercury_output and any( + self._point_uses_group_center(definition.get(key, [])) + for key in [ + "point_1_atoms", + "point_2_atoms", + "point_3_atoms", + "point_4_atoms", + ] + ): + mercury_torsion_values[label] = self.point_torsion( + coords, + definition.get("point_1_atoms", []), + definition.get("point_2_atoms", []), + definition.get("point_3_atoms", []), + definition.get("point_4_atoms", []), + definition.get("point_1_symmetry") or None, + definition.get("point_2_symmetry") or None, + definition.get("point_3_symmetry") or None, + definition.get("point_4_symmetry") or None, + group_center_round_dp=3, + ) + + for index, definition in enumerate(plane_distance_definitions or []): + if not isinstance(definition, dict): + continue + label = definition.get("label", f"Point_Plane_Distance_{index + 1}") + plane_distance_values[label] = self.point_plane_distance( + coords, + definition.get("point_atoms", []), + definition.get("plane_atoms", []), + definition.get("point_symmetry") or None, + definition.get("plane_symmetry") or None, + ) + if mercury_output and self._point_uses_group_center( + definition.get("point_atoms", []) + ): + mercury_plane_distance_values[label] = self.point_plane_distance( + coords, + definition.get("point_atoms", []), + definition.get("plane_atoms", []), + definition.get("point_symmetry") or None, + definition.get("plane_symmetry") or None, + group_center_round_dp=3, + ) + + self._append_measurements_csv( + "point_geometry_distances.csv", + structure_number, + distance_values, + results_path, + ) + self._append_measurements_csv( + "point_geometry_angles.csv", + structure_number, + angle_values, + results_path, + ) + self._append_measurements_csv( + "point_geometry_torsions.csv", + structure_number, + torsion_values, + results_path, + ) + self._append_measurements_csv( + "point_geometry_plane_distances.csv", + structure_number, + plane_distance_values, + results_path, + ) + + if mercury_output: + self._append_measurements_csv( + "point_geometry_distances_mercury.csv", + structure_number, + mercury_distance_values, + results_path, + ) + self._append_measurements_csv( + "point_geometry_angles_mercury.csv", + structure_number, + mercury_angle_values, + results_path, + ) + self._append_measurements_csv( + "point_geometry_torsions_mercury.csv", + structure_number, + mercury_torsion_values, + results_path, + ) + self._append_measurements_csv( + "point_geometry_plane_distances_mercury.csv", + structure_number, + mercury_plane_distance_values, + results_path, + ) + + def point_group_distance( + self, + coords: dict, + atom_list_1: list, + atom_list_2: list, + symmetry_1: str = None, + symmetry_2: str = None, + ) -> float: + """Calculates the distance in Angstroms between the centroids of two atom groups. + + Centroids are calculated in fractional coordinates then converted to Cartesian + space using the unit cell orthonormalisation matrix before computing the + Euclidean distance. + + Args: + coords (dict): fractional coordinates from LST_Read.extract_atom_coordinates() + atom_list_1 (list): atom labels for the first group + atom_list_2 (list): atom labels for the second group + symmetry_1 (str): optional symmetry operation string for centroid 1 + symmetry_2 (str): optional symmetry operation string for centroid 2 + + Returns: + distance (float): distance in Angstroms between the two centroids, + or 0.0 if either centroid could not be calculated + """ + + return self.point_distance( + coords, atom_list_1, atom_list_2, symmetry_1, symmetry_2 + ) + + def find_point_group_distance( + self, + file_name: str, + atom_list_1: list, + atom_list_2: list, + symmetry_1: str = None, + symmetry_2: str = None, + ) -> float: + """Reads a .lst file and calculates the centroid-to-centroid distance. + + Args: + file_name (str): full path to the .lst file + atom_list_1 (list): atom labels for the first group + atom_list_2 (list): atom labels for the second group + symmetry_1 (str): optional symmetry operation string for centroid 1 + symmetry_2 (str): optional symmetry operation string for centroid 2 + + Returns: + distance (float): distance in Angstroms between the two centroids + """ + + data = self.lst_reader.read(file_name) + coords = self.lst_reader.extract_atom_coordinates(data) + return self.point_group_distance( + coords, atom_list_1, atom_list_2, symmetry_1, symmetry_2 + ) + + def analyse_point_group_distance( + self, + lst_name: str, + structure_number: int, + results_path: str, + atom_list_1: list, + atom_list_2: list, + label: str = "Centroid Distance", + symmetry_1: str = None, + symmetry_2: str = None, + ) -> None: + """Calculates the centroid distance between two atom groups and appends to a .csv. + + The output csv filename is derived from the label parameter. + + Args: + lst_name (str): full path to the .lst file + structure_number (int): structure number as the independent variable + results_path (str): full path to the output results directory + atom_list_1 (list): atom labels for the first group + atom_list_2 (list): atom labels for the second group + label (str): column header for the distance in the output csv; + also used to derive the csv filename + symmetry_1 (str): optional symmetry operation string for centroid 1 + symmetry_2 (str): optional symmetry operation string for centroid 2 + """ + + if lst_name == "": + logging.info(__name__ + " : Refinement failed, no structure to analyse") + return + + self.grab_cell(pathlib.Path(lst_name)) + + if self.bad_flag: + return + + distance = self.find_point_group_distance( + pathlib.Path(lst_name), atom_list_1, atom_list_2, symmetry_1, symmetry_2 + ) + + df = pd.DataFrame({"Structure": [structure_number], label: [distance]}) + csv_name = label.lower().replace(" ", "_") + ".csv" + + os.chdir(results_path) + try: + old_data = pd.read_csv(csv_name) + except FileNotFoundError: + df.to_csv(csv_name, index=None) + else: + new_df = pd.concat([old_data, df]) + new_df.to_csv(csv_name, index=None) diff --git a/cx_asap/post_refinement_analysis/pipelines/points_pipeline.py b/cx_asap/post_refinement_analysis/pipelines/points_pipeline.py new file mode 100644 index 0000000..1480a08 --- /dev/null +++ b/cx_asap/post_refinement_analysis/pipelines/points_pipeline.py @@ -0,0 +1,280 @@ +#!/usr/bin/env python3 + +################################################################################################### +# --------------------------------------CX-ASAP: points_pipeline-----------------------------------# +# ---Authors: Amy J. Thompson, Kate M. Smith, Daniel J. Eriksson, Jack K. Clegg & Jason R. Price---# +# -----------------------------------Python Implementation by AJT----------------------------------# +# -----------------------------------Project Design by JRP and JKC---------------------------------# +# --------------------------------Valuable Coding Support by KMS & DJE-----------------------------# +################################################################################################### + +# ----------Required Modules----------# + +from system_files.utils import Config, Directory_Browse, Grapher +from post_refinement_analysis.modules.points import PointGeometryEngine +import os +import pathlib +import pandas as pd +import logging + +# ----------Class Definition----------# + + +class PointsPipeline: + def __init__(self) -> None: + """Initialises the class + + Sets up the yaml parameters input by the user + + Also defines the location for the system yaml file + + which stores a yaml of code-only parameters accessible throughout + + the software package + """ + + config = Config() + + self.cfg = config.cfg + self.sys = config.sys + self.conf_path = config.conf_path + self.sys_path = config.sys_path + + def create_numbered_results_directory( + self, base_directory: str, folder_name: str + ) -> pathlib.Path: + """Creates a numbered results directory inside base_directory/folder_name.""" + + base_path = pathlib.Path(base_directory) + results_root = base_path / folder_name + results_root.mkdir(exist_ok=True) + + existing_numbers = [] + for item in results_root.iterdir(): + if item.is_dir(): + try: + existing_numbers.append(int(item.name)) + except ValueError: + continue + + next_number = max(existing_numbers) + 1 if existing_numbers else 1 + results_path = results_root / str(next_number) + results_path.mkdir() + return results_path + + def point_group_distance_analysis( + self, + working_directory: str, + atom_list_1: list, + atom_list_2: list, + results_directory: str, + label: str = "Centroid Distance", + symmetry_1: str = None, + symmetry_2: str = None, + ) -> None: + """Calculates centroid-to-centroid distances for a series of .lst files + in separate folders within a common parent folder. + + Produces a scatter graph of centroid distance vs structure number. + + Args: + working_directory (str): full path to the parent folder containing + folders with .lst files + atom_list_1 (list): atom labels for the first point group + atom_list_2 (list): atom labels for the second point group + results_directory (str): full path to the output directory + label (str): column header for the output csv and graph y-axis + symmetry_1 (str): optional symmetry operation string for centroid 1 + symmetry_2 (str): optional symmetry operation string for centroid 2 + """ + + logging.info( + "Running centroid distance analysis with " + f"working_directory={working_directory}, " + f"results_directory={results_directory}, " + f"atom_list_1={atom_list_1}, " + f"atom_list_2={atom_list_2}, " + f"symmetry_1={symmetry_1}, " + f"symmetry_2={symmetry_2}" + ) + + centroid = PointGeometryEngine() + tree = Directory_Browse(working_directory) + + for index, item in enumerate(tree.directories): + tree.enter_directory(item, ".lst") + centroid.analyse_point_group_distance( + tree.item_file, + index + 1, + results_directory, + atom_list_1, + atom_list_2, + label, + symmetry_1, + symmetry_2, + ) + tree.exit_directory() + + os.chdir(results_directory) + csv_name = label.lower().replace(" ", "_") + ".csv" + + try: + full_data = pd.read_csv(csv_name) + x = full_data["Structure"] + y = full_data[label] + graph = Grapher() + graph.single_scatter_graph( + x, + y, + "Structure Number", + r"Distance ($\AA$)", + label, + csv_name.replace(".csv", ".png"), + ) + except FileNotFoundError: + logging.error(f"No {csv_name} file found...") + + def point_geometry_analysis( + self, + working_directory: str, + results_directory: str, + distance_definitions: "list[dict]" = None, + angle_definitions: "list[dict]" = None, + torsion_definitions: "list[dict]" = None, + plane_distance_definitions: "list[dict]" = None, + mercury_output: bool = False, + ) -> None: + """Calculates point-based geometry values across multiple .lst files. + + Supported metrics are: + - 2-point distances + - 3-point angles + - 4-point torsions + - point-to-plane distances + """ + + logging.info( + "Running point geometry analysis with " + f"working_directory={working_directory}, " + f"results_directory={results_directory}, " + f"distance_definitions={len(distance_definitions or [])}, " + f"angle_definitions={len(angle_definitions or [])}, " + f"torsion_definitions={len(torsion_definitions or [])}, " + f"plane_distance_definitions={len(plane_distance_definitions or [])}" + ) + + centroid = PointGeometryEngine() + tree = Directory_Browse(working_directory) + processed_structures = 0 + + for index, item in enumerate(tree.directories): + tree.enter_directory(item, ".lst") + centroid.analyse_point_geometry( + tree.item_file, + index + 1, + results_directory, + distance_definitions, + angle_definitions, + torsion_definitions, + plane_distance_definitions, + mercury_output=mercury_output, + ) + processed_structures += 1 + tree.exit_directory() + + os.chdir(results_directory) + + self._graph_point_geometry_csv( + "point_geometry_distances.csv", + "Point Geometry Distances", + r"Distance ($\AA$)", + "point_geometry_distances.png", + ) + if mercury_output: + self._graph_point_geometry_csv( + "point_geometry_distances_mercury.csv", + "Point Geometry Distances (Mercury)", + r"Distance ($\AA$)", + "point_geometry_distances_mercury.png", + ) + self._graph_point_geometry_csv( + "point_geometry_angles.csv", + "Point Geometry Angles", + "Angle($^\\circ$)", + "point_geometry_angles.png", + ) + if mercury_output: + self._graph_point_geometry_csv( + "point_geometry_angles_mercury.csv", + "Point Geometry Angles (Mercury)", + "Angle($^\\circ$)", + "point_geometry_angles_mercury.png", + ) + self._graph_point_geometry_csv( + "point_geometry_torsions.csv", + "Point Geometry Torsions", + "Angle($^\\circ$)", + "point_geometry_torsions.png", + ) + if mercury_output: + self._graph_point_geometry_csv( + "point_geometry_torsions_mercury.csv", + "Point Geometry Torsions (Mercury)", + "Angle($^\\circ$)", + "point_geometry_torsions_mercury.png", + ) + self._graph_point_geometry_csv( + "point_geometry_plane_distances.csv", + "Point to Plane Distances", + r"Distance ($\AA$)", + "point_geometry_plane_distances.png", + ) + if mercury_output: + self._graph_point_geometry_csv( + "point_geometry_plane_distances_mercury.csv", + "Point to Plane Distances (Mercury)", + r"Distance ($\AA$)", + "point_geometry_plane_distances_mercury.png", + ) + + logging.info( + "Point geometry analysis finished with " + f"processed_structures={processed_structures}, " + f"distances_csv={os.path.exists('point_geometry_distances.csv')}, " + f"angles_csv={os.path.exists('point_geometry_angles.csv')}, " + f"torsions_csv={os.path.exists('point_geometry_torsions.csv')}, " + "plane_distances_csv=" + f"{os.path.exists('point_geometry_plane_distances.csv')}, " + "mercury_output=" + f"{mercury_output}" + ) + + def _graph_point_geometry_csv( + self, csv_name: str, graph_title: str, y_axis_title: str, figure_name: str + ) -> None: + """Creates a summary scatter graph for one point-geometry CSV file.""" + + try: + full_data = pd.read_csv(csv_name) + except FileNotFoundError: + return + + if "Structure" not in full_data.columns: + return + + y_cols = [col for col in full_data.columns if col != "Structure"] + if not y_cols: + return + + x = full_data["Structure"] + y_data = [list(full_data[col]) for col in y_cols] + graph = Grapher() + graph.single_scatter_graph( + x, + y_data, + "Structure Number", + y_axis_title, + graph_title, + figure_name, + y_series_title=y_cols if len(y_cols) > 1 else None, + ) diff --git a/cx_asap/system_files/parameter.yaml b/cx_asap/system_files/parameter.yaml index bcbb452..09c76cd 100644 --- a/cx_asap/system_files/parameter.yaml +++ b/cx_asap/system_files/parameter.yaml @@ -67,11 +67,43 @@ pipeline-variable-position: - minimum_fraction_of_indexed_spots - mapping_step_size - total_angle - - calculate_centroid_distance - - centroid_1_atoms - - centroid_2_atoms - - centroid_1_symmetry - - centroid_2_symmetry + - calculate_point_group_distance + - point_group_1_atoms + - point_group_2_atoms + - point_group_1_symmetry + - point_group_2_symmetry + +module-point-geometry: + - lst_file_location + - point_geometry_distances + - point_geometry_angles + - point_geometry_torsions + - point_geometry_plane_distances + - mercury_output + +pipeline-point-geometry: + - experiment_location + - point_geometry_distances + - point_geometry_angles + - point_geometry_torsions + - point_geometry_plane_distances + - mercury_output + - structural_analysis_hbonds + - ADP_analysis + - experiment_location + - reference_unit_cell + - atoms_for_plane + - wedge_angles + - min_pixels + - spot_maximum_centroid + - signal_pixel + - sepmin + - mapping_step_size + - calculate_point_group_distance + - point_group_1_atoms + - point_group_2_atoms + - point_group_1_symmetry + - point_group_2_symmetry pipeline-general: - experiment_location @@ -313,13 +345,6 @@ module-rotation-planes: - lst_file_location - calculate_interplane_angle -module-centroids: - - lst_file_location - - centroid_1_atoms - - centroid_2_atoms - - centroid_1_symmetry - - centroid_2_symmetry - module-point-geometry: - lst_file_location - point_geometry_distances @@ -350,13 +375,6 @@ module-cif-analysis: - calculate_interplane_angle - lst_file_location -pipeline-centroids: - - experiment_location - - centroid_1_atoms - - centroid_2_atoms - - centroid_1_symmetry - - centroid_2_symmetry - pipeline-point-geometry: - experiment_location - point_geometry_distances @@ -430,11 +448,11 @@ pipeline-position-analysis: - signal_pixel - sepmin - mapping_step_size - - calculate_centroid_distance - - centroid_1_atoms - - centroid_2_atoms - - centroid_1_symmetry - - centroid_2_symmetry + - calculate_point_group_distance + - point_group_1_atoms + - point_group_2_atoms + - point_group_1_symmetry + - point_group_2_symmetry pipeline-temperature-analysis: - cif_parameters diff --git a/tests/test_yamls.py b/tests/test_yamls.py index 7bf8c54..3d76393 100644 --- a/tests/test_yamls.py +++ b/tests/test_yamls.py @@ -72,11 +72,51 @@ def setUp(self): "minimum_fraction_of_indexed_spots", "mapping_step_size", "total_angle", - "calculate_centroid_distance", - "centroid_1_atoms", - "centroid_2_atoms", - "centroid_1_symmetry", - "centroid_2_symmetry", + "calculate_point_group_distance", + "point_group_1_atoms", + "point_group_2_atoms", + "point_group_1_symmetry", + "point_group_2_symmetry", + ], + "module-point-geometry": [ + "lst_file_location", + "point_geometry_distances", + "point_geometry_angles", + "point_geometry_torsions", + "point_geometry_plane_distances", + "mercury_output", + ], + "pipeline-point-geometry": [ + "experiment_location", + "point_geometry_distances", + "point_geometry_angles", + "point_geometry_torsions", + "point_geometry_plane_distances", + "mercury_output", + ], + "pipeline-position-analysis": [ + "cif_parameters", + "atoms_for_analysis", + "reference_plane", + "structural_analysis_bonds", + "structural_analysis_angles", + "structural_analysis_torsions", + "structural_analysis_hbonds", + "ADP_analysis", + "experiment_location", + "reference_unit_cell", + "atoms_for_plane", + "wedge_angles", + "min_pixels", + "spot_maximum_centroid", + "signal_pixel", + "sepmin", + "mapping_step_size", + "calculate_point_group_distance", + "point_group_1_atoms", + "point_group_2_atoms", + "point_group_1_symmetry", + "point_group_2_symmetry", ], "pipeline-general": [ "experiment_location", @@ -347,11 +387,11 @@ def setUp(self): "signal_pixel", "sepmin", "mapping_step_size", - "calculate_centroid_distance", - "centroid_1_atoms", - "centroid_2_atoms", - "centroid_1_symmetry", - "centroid_2_symmetry", + "calculate_point_group_distance", + "point_group_1_atoms", + "point_group_2_atoms", + "point_group_1_symmetry", + "point_group_2_symmetry", ], "pipeline-temperature-analysis": [ "cif_parameters", @@ -443,11 +483,51 @@ def setUp(self): "minimum_fraction_of_indexed_spots", "mapping_step_size", "total_angle", - "calculate_centroid_distance", - "centroid_1_atoms", - "centroid_2_atoms", - "centroid_1_symmetry", - "centroid_2_symmetry", + "calculate_point_group_distance", + "point_group_1_atoms", + "point_group_2_atoms", + "point_group_1_symmetry", + "point_group_2_symmetry", + ], + [ + "lst_file_location", + "point_geometry_distances", + "point_geometry_angles", + "point_geometry_torsions", + "point_geometry_plane_distances", + "mercury_output", + ], + [ + "experiment_location", + "point_geometry_distances", + "point_geometry_angles", + "point_geometry_torsions", + "point_geometry_plane_distances", + "mercury_output", + ], + [ + "cif_parameters", + "atoms_for_analysis", + "reference_plane", + "structural_analysis_bonds", + "structural_analysis_angles", + "structural_analysis_torsions", + "structural_analysis_hbonds", + "ADP_analysis", + "experiment_location", + "reference_unit_cell", + "atoms_for_plane", + "wedge_angles", + "min_pixels", + "spot_maximum_centroid", + "signal_pixel", + "sepmin", + "mapping_step_size", + "calculate_point_group_distance", + "point_group_1_atoms", + "point_group_2_atoms", + "point_group_1_symmetry", + "point_group_2_symmetry", ], [ "experiment_location", @@ -684,30 +764,6 @@ def setUp(self): "experiment_location", "reference_unit_cell", ], - [ - "cif_parameters", - "atoms_for_analysis", - "reference_plane", - "structural_analysis_bonds", - "structural_analysis_angles", - "structural_analysis_torsions", - "structural_analysis_hbonds", - "ADP_analysis", - "experiment_location", - "reference_unit_cell", - "atoms_for_plane", - "wedge_angles", - "min_pixels", - "spot_maximum_centroid", - "signal_pixel", - "sepmin", - "mapping_step_size", - "calculate_centroid_distance", - "centroid_1_atoms", - "centroid_2_atoms", - "centroid_1_symmetry", - "centroid_2_symmetry", - ], [ "cif_parameters", "atoms_for_analysis",