diff --git a/environment.yml b/environment.yml index 247f36f1c..9da35a67c 100644 --- a/environment.yml +++ b/environment.yml @@ -12,7 +12,6 @@ dependencies: - lomap2>=3.2.1 - networkx - numpy - - openfe-analysis>=0.4.0 # min pin https://github.com/OpenFreeEnergy/openfe/issues/1834#issuecomment-3920079481, no max to check issues with new versions - openff-interchange-base >=0.5.0,!= 0.5.1 # https://github.com/openforcefield/openff-interchange/issues/1450 and https://github.com/OpenFreeEnergy/openfe/pull/1901 - openff-nagl-base >=0.3.3 - openff-nagl-models>=0.1.2 @@ -54,6 +53,7 @@ dependencies: - threadpoolctl - pip: - git+https://github.com/OpenFreeEnergy/gufe@main + - git+https://github.com/OpenFreeEnergy/openfe_analysis@main - run_constrained: # drop this pin when handled upstream in espaloma-feedstock - smirnoff99frosst>=1.1.0.1 #https://github.com/openforcefield/smirnoff99Frosst/issues/109 diff --git a/src/openfe/protocols/openmm_septop/base_units.py b/src/openfe/protocols/openmm_septop/base_units.py index 821c1931a..7d1724013 100644 --- a/src/openfe/protocols/openmm_septop/base_units.py +++ b/src/openfe/protocols/openmm_septop/base_units.py @@ -20,6 +20,9 @@ from typing import Any, Literal, Optional import gufe +import matplotlib.pyplot as plt +import MDAnalysis as mda +import numpy as np import numpy.typing as npt import openmm import openmmtools @@ -43,6 +46,7 @@ ThermodynamicState, create_thermodynamic_state_protocol, ) +from rdkit import Chem import openfe from openfe.protocols.openmm_afe.equil_afe_settings import ( @@ -1503,11 +1507,302 @@ def _analyze_multistate_energies( reporter.close() return analyzer.unit_results_dict + @classmethod + def _run_complex_analysis( + cls, + ds, + pdb_file: pathlib.Path, + skip: int, + ligand_A_indices: list[int], + ligand_B_indices: list[int], + rdmol_A: Chem.Mol, + rdmol_B: Chem.Mol, + protein_selection: str = "protein and name CA", + ) -> dict[str, Any]: + """ + Run structural analysis for the complex phase. + + Parameters + ---------- + ds : netCDF4.Dataset + Open NetCDF dataset for the multistate trajectory. + pdb_file : pathlib.Path + Path to the subsampled PDB file. + skip : int + Frame stride for analysis. + ligand_A_indices : list[int] + Atom indices of ligand A in the subsampled system. + ligand_B_indices : list[int] + Atom indices of ligand B in the subsampled system. + rdmol_A : Chem.Mol + RDKit molecule for ligand A, used for symmetry-corrected RMSD. + rdmol_B : Chem.Mol + RDKit molecule for ligand B, used for symmetry-corrected RMSD. + protein_selection : str + MDAnalysis selection string for the protein atoms used for + alignment and RMSD calculations. Default: "protein and name CA". + + Returns + ------- + dict[str, Any] + Per-state lists for ligand RMSD, COM drift, protein 2D RMSD, + and a single time_ps array. + """ + from openfe_analysis.rmsd import ( + LigandCOMDrift, + Protein2DRMSD, + SymmetryCorrectedLigandRMSD, + ) + from openfe_analysis.utils.apply_transformations import ( + apply_complex_alignment_transformations, + ) + from openfe_analysis.utils.universe_utils import create_universe_single_state + + n_lambda = ds.dimensions["state"].size + data: dict[str, Any] = { + "ligand_A_RMSD": [], + "ligand_B_RMSD": [], + "ligand_A_COM_drift": [], + "ligand_B_COM_drift": [], + "protein_2D_RMSD": [], + "time_ps": None, + } + u_top = mda.Universe(pdb_file) + for state_idx in range(n_lambda): + universe = create_universe_single_state(u_top._topology, ds, state=state_idx) + prot = universe.select_atoms(protein_selection) + lig_A = universe.atoms[ligand_A_indices] + lig_B = universe.atoms[ligand_B_indices] + apply_complex_alignment_transformations(universe, protein=prot, ligands=[lig_A, lig_B]) + + if prot: + prot_rmsd2d = Protein2DRMSD(prot).run(step=skip) + data["protein_2D_RMSD"].append(prot_rmsd2d.results.rmsd2d) + + for label, lig, rdmol in [ + ("ligand_A", lig_A, rdmol_A), + ("ligand_B", lig_B, rdmol_B), + ]: + lig_rmsd = SymmetryCorrectedLigandRMSD(lig, rdmol=rdmol).run(step=skip) + data[f"{label}_RMSD"].append(lig_rmsd.results.rmsd) + + lig_drift = LigandCOMDrift(lig).run(step=skip) + data[f"{label}_COM_drift"].append(lig_drift.results.com_drift) + + if data["time_ps"] is None: + data["time_ps"] = ( + np.arange(len(universe.trajectory))[::skip] * universe.trajectory.dt + ) + + return data + + @classmethod + def _run_solvent_analysis( + cls, + ds, + pdb_file: pathlib.Path, + skip: int, + ligand_A_indices: list[int], + ligand_B_indices: list[int], + rdmol_A: Chem.Mol, + rdmol_B: Chem.Mol, + ) -> dict[str, Any]: + """ + Run structural analysis for the solvent phase. + + Parameters + ---------- + ds : netCDF4.Dataset + Open NetCDF dataset for the multistate trajectory. + pdb_file : pathlib.Path + Path to the subsampled PDB file. + skip : int + Frame stride for analysis. + ligand_A_indices : list[int] + Atom indices of ligand A in the subsampled system. + ligand_B_indices : list[int] + Atom indices of ligand B in the subsampled system. + rdmol_A : Chem.Mol + RDKit molecule for ligand A, used for symmetry-corrected RMSD. + rdmol_B : Chem.Mol + RDKit molecule for ligand B, used for symmetry-corrected RMSD. + + Returns + ------- + dict[str, Any] + Per-state lists for ligand RMSD and a single + time_ps array. + """ + from openfe_analysis.rmsd import SymmetryCorrectedLigandRMSD + from openfe_analysis.utils.apply_transformations import ( + apply_ligand_alignment_transformations, + ) + from openfe_analysis.utils.universe_utils import create_universe_single_state + + n_lambda = ds.dimensions["state"].size + data: dict[str, Any] = { + "ligand_A_RMSD": [], + "ligand_B_RMSD": [], + "time_ps": None, + } + u_top = mda.Universe(pdb_file) + for state_idx in range(n_lambda): + for label, indices, rdmol in [ + ("ligand_A", ligand_A_indices, rdmol_A), + ("ligand_B", ligand_B_indices, rdmol_B), + ]: + universe = create_universe_single_state(u_top._topology, ds, state=state_idx) + lig = universe.atoms[indices] + apply_ligand_alignment_transformations(universe, ligand=lig) + + lig_rmsd = SymmetryCorrectedLigandRMSD(lig, rdmol=rdmol).run(step=skip) + data[f"{label}_RMSD"].append(lig_rmsd.results.rmsd) + + if data["time_ps"] is None: + data["time_ps"] = ( + np.arange(len(universe.trajectory))[::skip] * universe.trajectory.dt + ) + + return data + + @classmethod + def _structural_analysis( + cls, + pdb_file: pathlib.Path, + trj_file: pathlib.Path, + output_directory: pathlib.Path, + dry: bool, + simtype: str, + ligand_A_indices: list[int], + ligand_B_indices: list[int], + rdmol_A: Chem.Mol, + rdmol_B: Chem.Mol, + protein_selection: str = "protein and name CA", + ) -> dict[str, str | pathlib.Path]: + """ + Run structural analysis using ``openfe-analysis``. + + Parameters + ---------- + pdb_file : pathlib.Path + Path to the PDB file. + trj_file : pathlib.Path + Path to the trajectory NetCDF file. + output_directory : pathlib.Path + The output directory where plots and the data NPZ file + will be stored. + dry : bool + Whether or not we are running a dry run. + simtype : str + Either ``"complex"`` or ``"solvent"``. Controls whether protein + analyses are run and how alignment is applied. + ligand_A_indices : list[int] + Atom indices of ligand A in the subsampled system. + ligand_B_indices : list[int] + Atom indices of ligand B in the subsampled system. + rdmol_A : Chem.Mol + RDKit molecule for ligand A, used for symmetry-corrected RMSD. + rdmol_B : Chem.Mol + RDKit molecule for ligand B, used for symmetry-corrected RMSD. + protein_selection : str + MDAnalysis selection string for the protein atoms used for + alignment and RMSD calculations in the complex phase. + Ignored for the solvent phase. Default "protein and name CA". + + Returns + ------- + dict[str, str | pathlib.Path] + Dictionary containing either the path to the NPZ file with the + structural data, or the analysis error. + """ + import netCDF4 as nc + from openfe_analysis.utils import plotting + + try: + with nc.Dataset(trj_file) as ds: + if hasattr(ds, "PositionInterval"): + n_frames = len(range(0, ds.dimensions["iteration"].size, ds.PositionInterval)) + else: + n_frames = ds.dimensions["iteration"].size + skip = max(n_frames // 500, 1) + + if simtype == "complex": + data = cls._run_complex_analysis( + ds, + pdb_file, + skip, + ligand_A_indices, + ligand_B_indices, + rdmol_A, + rdmol_B, + protein_selection=protein_selection, + ) + else: + data = cls._run_solvent_analysis( + ds, + pdb_file, + skip, + ligand_A_indices, + ligand_B_indices, + rdmol_A, + rdmol_B, + ) + + except Exception as e: + return {"structural_analysis_error": str(e)} + + # Generate relevant plots if not a dry run + if not dry: + time_ps = data.get("time_ps", []) + + if data.get("protein_2D_RMSD"): + fig = plotting.plot_2D_rmsd(data["protein_2D_RMSD"]) + fig.savefig(output_directory / "protein_2D_RMSD.png") + plt.close(fig) + + for label in ["ligand_A", "ligand_B"]: + if data.get(f"{label}_RMSD"): + fig = plotting.plot_ligand_RMSD(time_ps, data[f"{label}_RMSD"]) + fig.savefig(output_directory / f"{label}_RMSD.png") + plt.close(fig) + + if data.get(f"{label}_COM_drift"): + fig = plotting.plot_ligand_COM_drift(time_ps, data[f"{label}_COM_drift"]) + fig.savefig(output_directory / f"{label}_COM_drift.png") + plt.close(fig) + + # Write out an NPZ with all the relevant analysis data + npz_file = output_directory / "structural_analysis.npz" + if simtype == "complex": + np.savez_compressed( + npz_file, + ligand_A_RMSD=np.asarray(data["ligand_A_RMSD"], dtype=np.float32), + ligand_B_RMSD=np.asarray(data["ligand_B_RMSD"], dtype=np.float32), + ligand_A_COM_drift=np.asarray(data["ligand_A_COM_drift"], dtype=np.float32), + ligand_B_COM_drift=np.asarray(data["ligand_B_COM_drift"], dtype=np.float32), + protein_2D_RMSD=np.asarray(data["protein_2D_RMSD"], dtype=np.float32), + time_ps=np.asarray(data["time_ps"], dtype=np.float32), + ) + else: + np.savez_compressed( + npz_file, + ligand_A_RMSD=np.asarray(data["ligand_A_RMSD"], dtype=np.float32), + ligand_B_RMSD=np.asarray(data["ligand_B_RMSD"], dtype=np.float32), + time_ps=np.asarray(data["time_ps"], dtype=np.float32), + ) + return {"structural_analysis": npz_file} + def run( self, *, trajectory: pathlib.Path, checkpoint: pathlib.Path, + pdb_file: pathlib.Path, + ligand_A_indices: list[int], + ligand_B_indices: list[int], + rdmol_A: Chem.Mol, + rdmol_B: Chem.Mol, + protein_selection: str = "protein and name CA", dry: bool = False, verbose: bool = True, scratch_basepath: pathlib.Path | None = None, @@ -1521,6 +1816,20 @@ def run( Path to the MultiStateReporter generated NetCDF file. checkpoint : pathlib.Path Path to the checkpoint file generated by MultiStateReporter. + pdb_file : pathlib.Path + Path to the subsampled PDB file. + ligand_A_indices : list[int] + Atom indices of ligand A in the subsampled system. + ligand_B_indices : list[int] + Atom indices of ligand B in the subsampled system. + rdmol_A : Chem.Mol + RDKit molecule for ligand A. + rdmol_B : Chem.Mol + RDKit molecule for ligand B. + protein_selection : str + MDAnalysis selection string for the protein atoms used for + alignment and RMSD calculations in the complex phase. + Ignored for the solvent phase. Default "protein and name CA". dry : bool Do a dry run of the calculation, creating all necessary hybrid system components (topology, system, sampler, etc...) but without @@ -1549,7 +1858,7 @@ def run( settings = self._get_settings() # Energies analysis - if verbose: + if self.verbose: self.logger.info("Analyzing energies") energy_analysis = self._analyze_multistate_energies( @@ -1560,7 +1869,23 @@ def run( dry=dry, ) - return energy_analysis + if self.verbose: + self.logger.info("Analyzing structural outputs") + + structural_analysis = self._structural_analysis( + pdb_file=pdb_file, + trj_file=trajectory, + output_directory=self.shared_basepath, + dry=dry, + simtype=self.simtype, + ligand_A_indices=ligand_A_indices, + ligand_B_indices=ligand_B_indices, + rdmol_A=rdmol_A, + rdmol_B=rdmol_B, + protein_selection=protein_selection, + ) + + return energy_analysis | structural_analysis def _execute( self, @@ -1579,14 +1904,31 @@ def _execute( trajectory = simulation.outputs["trajectory"] checkpoint = simulation.outputs["checkpoint"] + alchem_comps = self._inputs["alchemical_components"] + rdmol_A = alchem_comps["stateA"][0].to_rdkit() + rdmol_B = alchem_comps["stateB"][0].to_rdkit() + + # Remap ligand indices from full system to subsampled system + selection_indices = np.array(setup.outputs["selection_indices"]) + ligand_A_indices = np.where(np.isin(selection_indices, setup.outputs["ligand_A_indices"]))[ + 0 + ].tolist() + ligand_B_indices = np.where(np.isin(selection_indices, setup.outputs["ligand_B_indices"]))[ + 0 + ].tolist() + outputs = self.run( trajectory=trajectory, checkpoint=checkpoint, + pdb_file=setup.outputs["subsampled_pdb_structure"], + ligand_A_indices=ligand_A_indices, + ligand_B_indices=ligand_B_indices, + rdmol_A=rdmol_A, + rdmol_B=rdmol_B, scratch_basepath=ctx.scratch, shared_basepath=ctx.shared, ) - # We re-include things here to make life easier when gathering results if self.simtype == "complex": previous_outputs = { "standard_state_correction_A": setup.outputs["standard_state_correction_A"], diff --git a/src/openfe/protocols/openmm_septop/equil_septop_method.py b/src/openfe/protocols/openmm_septop/equil_septop_method.py index 05d84a1bf..2f3fd8443 100644 --- a/src/openfe/protocols/openmm_septop/equil_septop_method.py +++ b/src/openfe/protocols/openmm_septop/equil_septop_method.py @@ -577,6 +577,9 @@ def _create( analysis = unit_classes[phase]["analysis"]( protocol=self, + stateA=stateA, + stateB=stateB, + alchemical_components=alchem_comps, setup=setup, simulation=simulation, generation=0, diff --git a/src/openfe/protocols/openmm_septop/septop_units.py b/src/openfe/protocols/openmm_septop/septop_units.py index 6ad261837..12fde00e5 100644 --- a/src/openfe/protocols/openmm_septop/septop_units.py +++ b/src/openfe/protocols/openmm_septop/septop_units.py @@ -855,6 +855,8 @@ def run( "restraint_geometry_B": restraint_geom_B.model_dump(), "selection_indices": selection_indices, "subsampled_pdb_structure": sub_pdb_structure, + "ligand_A_indices": atom_indices_AB_A, + "ligand_B_indices": atom_indices_AB_B, } else: return { @@ -870,6 +872,8 @@ def run( "positions": equil_positions_AB, "selection_indices": selection_indices, "subsampled_pdb_structure": sub_pdb_structure, + "ligand_A_indices": atom_indices_AB_A, + "ligand_B_indices": atom_indices_AB_B, } @@ -1133,6 +1137,8 @@ def run( "standard_state_correction": corr.to("kilocalorie_per_mole"), "selection_indices": selection_indices, "subsampled_pdb_structure": sub_pdb_structure, + "ligand_A_indices": atom_indices_AB_A, + "ligand_B_indices": atom_indices_AB_B, } else: return { @@ -1146,6 +1152,8 @@ def run( "positions": positions_AB, "selection_indices": selection_indices, "subsampled_pdb_structure": sub_pdb_structure, + "ligand_A_indices": atom_indices_AB_A, + "ligand_B_indices": atom_indices_AB_B, } diff --git a/src/openfe/tests/protocols/openmm_septop/test_septop_protocol_analysis.py b/src/openfe/tests/protocols/openmm_septop/test_septop_protocol_analysis.py new file mode 100644 index 000000000..5d57acf4e --- /dev/null +++ b/src/openfe/tests/protocols/openmm_septop/test_septop_protocol_analysis.py @@ -0,0 +1,214 @@ +import pathlib + +import numpy as np +import pooch +import pytest +from rdkit.Chem import SDMolSupplier + +from openfe.data._registry import POOCH_CACHE +from openfe.protocols.openmm_septop.base_units import BaseSepTopAnalysisUnit + +pooch_septop_structural = pooch.create( + path=POOCH_CACHE, + base_url="https://zenodo.org/records/20507106/files/", + registry={ + "septop_structural_results.zip": "md5:cffc193dacb5ad8d26c5467ae05b1a00", + }, +) + + +@pytest.fixture(scope="session") +def septop_structural_results_dir(): + pooch_septop_structural.fetch("septop_structural_results.zip", processor=pooch.Unzip()) + return pathlib.Path( + POOCH_CACHE / "septop_structural_results.zip.unzip/septop_structural_results" + ) + + +@pytest.fixture(scope="session") +def septop_complex_data(septop_structural_results_dir): + base = septop_structural_results_dir / "complex" + return { + "pdb": base / "alchemical_system.pdb", + "nc": base / "complex.nc", + "ligand_A_indices": np.load(base / "ligand_A_indices.npy").tolist(), + "ligand_B_indices": np.load(base / "ligand_B_indices.npy").tolist(), + "rdmol_A": SDMolSupplier(str(base / "ligand_A.sdf"), removeHs=False)[0], + "rdmol_B": SDMolSupplier(str(base / "ligand_B.sdf"), removeHs=False)[0], + } + + +@pytest.fixture(scope="session") +def septop_solvent_data(septop_structural_results_dir): + base = septop_structural_results_dir / "solvent" + return { + "pdb": base / "alchemical_system.pdb", + "nc": base / "solvent.nc", + "ligand_A_indices": np.load(base / "ligand_A_indices.npy").tolist(), + "ligand_B_indices": np.load(base / "ligand_B_indices.npy").tolist(), + "rdmol_A": SDMolSupplier(str(base / "ligand_A.sdf"), removeHs=False)[0], + "rdmol_B": SDMolSupplier(str(base / "ligand_B.sdf"), removeHs=False)[0], + } + + +@pytest.fixture(scope="session") +def complex_structural_analysis_result(septop_complex_data, tmp_path_factory): + d = septop_complex_data + tmp = tmp_path_factory.mktemp("complex_structural_analysis") + result = BaseSepTopAnalysisUnit._structural_analysis( + pdb_file=d["pdb"], + trj_file=d["nc"], + output_directory=tmp, + dry=False, + simtype="complex", + ligand_A_indices=d["ligand_A_indices"], + ligand_B_indices=d["ligand_B_indices"], + rdmol_A=d["rdmol_A"], + rdmol_B=d["rdmol_B"], + ) + return result, tmp + + +@pytest.fixture(scope="session") +def solvent_structural_analysis_result(septop_solvent_data, tmp_path_factory): + d = septop_solvent_data + tmp = tmp_path_factory.mktemp("solvent_structural_analysis") + result = BaseSepTopAnalysisUnit._structural_analysis( + pdb_file=d["pdb"], + trj_file=d["nc"], + output_directory=tmp, + dry=False, + simtype="solvent", + ligand_A_indices=d["ligand_A_indices"], + ligand_B_indices=d["ligand_B_indices"], + rdmol_A=d["rdmol_A"], + rdmol_B=d["rdmol_B"], + ) + return result, tmp + + +class TestComplexStructuralAnalysis: + def test_npz_written(self, complex_structural_analysis_result): + result, _ = complex_structural_analysis_result + assert "structural_analysis" in result + assert result["structural_analysis"].exists() + + def test_npz_keys_values_shape(self, complex_structural_analysis_result): + result, _ = complex_structural_analysis_result + npz = np.load(result["structural_analysis"]) + expected_keys = { + "ligand_A_RMSD", + "ligand_B_RMSD", + "ligand_A_COM_drift", + "ligand_B_COM_drift", + "protein_2D_RMSD", + "time_ps", + } + assert set(npz.files) == expected_keys + + # First frame RMSD should be zero (reference frame) + assert npz["ligand_A_RMSD"][0][0] == pytest.approx(0.0, abs=1e-5) + assert npz["ligand_B_RMSD"][0][0] == pytest.approx(0.0, abs=1e-5) + assert npz["ligand_A_COM_drift"][0][0] == pytest.approx(0.0, abs=1e-5) + assert npz["ligand_B_COM_drift"][0][0] == pytest.approx(0.0, abs=1e-5) + + # Time should start at zero + assert npz["time_ps"][0] == pytest.approx(0.0, abs=1e-5) + + # All RMSD and COM drift values should be non-negative + for state in npz["ligand_A_RMSD"]: + assert np.all(state >= 0) + for state in npz["ligand_B_RMSD"]: + assert np.all(state >= 0) + for state in npz["ligand_A_COM_drift"]: + assert np.all(state >= 0) + for state in npz["ligand_B_COM_drift"]: + assert np.all(state >= 0) + + # Should be 13 lambda windows + n_lambda = 13 + for key in expected_keys - {"time_ps"}: + assert len(npz[key]) == n_lambda + + def test_plots_written(self, complex_structural_analysis_result): + result, tmp = complex_structural_analysis_result + expected_plots = { + "ligand_A_RMSD.png", + "ligand_B_RMSD.png", + "ligand_A_COM_drift.png", + "ligand_B_COM_drift.png", + "protein_2D_RMSD.png", + } + written = {f.name for f in tmp.glob("*.png")} + assert written == expected_plots + + def test_bad_trajectory_returns_error_dict(self, septop_complex_data, tmp_path): + d = septop_complex_data + result = BaseSepTopAnalysisUnit._structural_analysis( + pdb_file=d["pdb"], + trj_file=tmp_path / "nonexistent.nc", + output_directory=tmp_path, + dry=True, + simtype="complex", + ligand_A_indices=d["ligand_A_indices"], + ligand_B_indices=d["ligand_B_indices"], + rdmol_A=d["rdmol_A"], + rdmol_B=d["rdmol_B"], + ) + + assert "structural_analysis_error" in result + assert "structural_analysis" not in result + + +class TestSolventStructuralAnalysis: + def test_npz_written(self, solvent_structural_analysis_result): + result, _ = solvent_structural_analysis_result + assert "structural_analysis" in result + assert result["structural_analysis"].exists() + + def test_npz_keys_values_shape(self, solvent_structural_analysis_result): + result, _ = solvent_structural_analysis_result + npz = np.load(result["structural_analysis"]) + expected_keys = {"ligand_A_RMSD", "ligand_B_RMSD", "time_ps"} + assert set(npz.files) == expected_keys + + # First frame RMSD should be zero (reference frame) + assert npz["ligand_A_RMSD"][0][0] == pytest.approx(0.0, abs=1e-5) + assert npz["ligand_B_RMSD"][0][0] == pytest.approx(0.0, abs=1e-5) + + # Time should start at zero + assert npz["time_ps"][0] == pytest.approx(0.0, abs=1e-5) + + # All RMSD and COM drift values should be non-negative + for state in npz["ligand_A_RMSD"]: + assert np.all(state >= 0) + for state in npz["ligand_B_RMSD"]: + assert np.all(state >= 0) + + # Should be 13 lambda windows + n_lambda = 13 + for key in expected_keys - {"time_ps"}: + assert len(npz[key]) == n_lambda + + def test_plots_written(self, solvent_structural_analysis_result): + result, tmp = solvent_structural_analysis_result + expected_plots = {"ligand_A_RMSD.png", "ligand_B_RMSD.png"} + written = {f.name for f in tmp.glob("*.png")} + assert written == expected_plots + + def test_bad_trajectory_returns_error_dict(self, septop_solvent_data, tmp_path): + d = septop_solvent_data + result = BaseSepTopAnalysisUnit._structural_analysis( + pdb_file=d["pdb"], + trj_file=tmp_path / "nonexistent.nc", + output_directory=tmp_path, + dry=True, + simtype="solvent", + ligand_A_indices=d["ligand_A_indices"], + ligand_B_indices=d["ligand_B_indices"], + rdmol_A=d["rdmol_A"], + rdmol_B=d["rdmol_B"], + ) + + assert "structural_analysis_error" in result + assert "structural_analysis" not in result diff --git a/src/openfe/tests/protocols/openmm_septop/test_septop_protocol_results.py b/src/openfe/tests/protocols/openmm_septop/test_septop_protocol_results.py index fded8adff..cfca396a6 100644 --- a/src/openfe/tests/protocols/openmm_septop/test_septop_protocol_results.py +++ b/src/openfe/tests/protocols/openmm_septop/test_septop_protocol_results.py @@ -76,6 +76,8 @@ def patcher(): ] ), "subsampled_pdb_structure": "subsampled.pdb", + "ligand_A_indices": [0], + "ligand_B_indices": [0], }, ), mock.patch( @@ -90,6 +92,8 @@ def patcher(): ] ), "subsampled_pdb_structure": "subsampled.pdb", + "ligand_A_indices": [0], + "ligand_B_indices": [0], }, ), mock.patch(