From 4b257b28a4e884f9a05fe24c14f5ffed95d1121a Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Thu, 18 Dec 2025 13:47:56 +0100 Subject: [PATCH 1/5] Fix stab at fixing multi chain RMSD analysis --- src/openfe_analysis/rmsd.py | 15 +++++++++------ src/openfe_analysis/transformations.py | 3 ++- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/src/openfe_analysis/rmsd.py b/src/openfe_analysis/rmsd.py index 7d2af13..2311784 100644 --- a/src/openfe_analysis/rmsd.py +++ b/src/openfe_analysis/rmsd.py @@ -7,6 +7,7 @@ import numpy as np import tqdm from MDAnalysis.analysis import rms +from MDAnalysis.transformations import make_whole, unwrap from numpy import typing as npt from .reader import FEReader @@ -42,15 +43,17 @@ def make_Universe(top: pathlib.Path, trj: nc.Dataset, state: int) -> mda.Univers if prot: # if there's a protein in the system: - # - make the protein not jump periodic images between frames + # - make the protein whole across periodic images between frames # - put the ligand in the closest periodic image as the protein # - align everything to minimise protein RMSD - nope = NoJump(prot) + make_whole_tr = make_whole(prot, compound="segments") + unwrap_tr = unwrap(prot) minnie = Minimiser(prot, ligand) align = Aligner(prot) u.trajectory.add_transformations( - nope, + make_whole_tr, + unwrap_tr, minnie, align, ) @@ -128,9 +131,9 @@ def gather_rms_data( # TODO: Some smart guard to avoid allocating a silly amount of memory? prot2d = np.empty((len(u.trajectory[::skip]), len(prot), 3), dtype=np.float32) - prot_start = prot.positions - # prot_weights = prot.masses / np.mean(prot.masses) - ligand_start = ligand.positions + # Would this copy be safer? + prot_start = prot.positions.copy() + ligand_start = ligand.positions.copy() ligand_initial_com = ligand.center_of_mass() ligand_weights = ligand.masses / np.mean(ligand.masses) diff --git a/src/openfe_analysis/transformations.py b/src/openfe_analysis/transformations.py index 1a4ef34..61db292 100644 --- a/src/openfe_analysis/transformations.py +++ b/src/openfe_analysis/transformations.py @@ -87,7 +87,8 @@ class Aligner(TransformationBase): def __init__(self, ref_ag: mda.AtomGroup): super().__init__() self.ref_idx = ref_ag.ix - self.ref_pos = ref_ag.positions + # Would this copy be safer? + self.ref_pos = ref_ag.positions.copy() self.weights = np.asarray(ref_ag.masses, dtype=np.float64) self.weights /= np.mean(self.weights) # normalise weights # remove COM shift from reference positions From 236ff7215383928351c8343f0f5afff08928e8f1 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Thu, 18 Dec 2025 14:56:54 +0100 Subject: [PATCH 2/5] Some updates --- src/openfe_analysis/cli.py | 46 ++++++++++++++++++++++++++----------- src/openfe_analysis/rmsd.py | 32 +++++++++++++++++++++++--- 2 files changed, 61 insertions(+), 17 deletions(-) diff --git a/src/openfe_analysis/cli.py b/src/openfe_analysis/cli.py index c45c540..437aad7 100644 --- a/src/openfe_analysis/cli.py +++ b/src/openfe_analysis/cli.py @@ -12,18 +12,36 @@ def cli(): @cli.command(name="RFE_analysis") -@click.argument( - "loc", - type=click.Path( - exists=True, readable=True, file_okay=False, dir_okay=True, path_type=pathlib.Path - ), +@click.option( + "--pdb", + type=click.Path(exists=True, readable=True, dir_okay=False, path_type=pathlib.Path), + required=True, + help="Path to the topology PDB file.", ) -@click.argument("output", type=click.Path(writable=True, dir_okay=False, path_type=pathlib.Path)) -def rfe_analysis(loc, output): - pdb = loc / "hybrid_system.pdb" - trj = loc / "simulation.nc" - - data = rmsd.gather_rms_data(pdb, trj) - - with click.open_file(output, "w") as f: - f.write(json.dumps(data)) +@click.option( + "--nc", + type=click.Path(exists=True, readable=True, dir_okay=False, path_type=pathlib.Path), + required=True, + help="Path to the NetCDF trajectory file.", +) +@click.option( + "--output", + type=click.Path(writable=True, dir_okay=False, path_type=pathlib.Path), + required=True, + help="Path to save the JSON results.", +) +def rfe_analysis(pdb: pathlib.Path, nc: pathlib.Path, output: pathlib.Path): + """ + Perform RMSD analysis for an RBFE simulation. + + Arguments: + pdb: path to the topology PDB file. + nc: path to the trajectory file (NetCDF format). + output: path to save the JSON results. + """ + # Run RMSD analysis + data = rmsd.gather_rms_data(pdb, nc) + + # Write results + with output.open("w") as f: + json.dump(data, f, indent=2) \ No newline at end of file diff --git a/src/openfe_analysis/rmsd.py b/src/openfe_analysis/rmsd.py index 2311784..1d9f29e 100644 --- a/src/openfe_analysis/rmsd.py +++ b/src/openfe_analysis/rmsd.py @@ -7,12 +7,31 @@ import numpy as np import tqdm from MDAnalysis.analysis import rms -from MDAnalysis.transformations import make_whole, unwrap +from MDAnalysis.transformations import unwrap +from MDAnalysis.lib.mdamath import make_whole from numpy import typing as npt from .reader import FEReader from .transformations import Aligner, Minimiser, NoJump +from MDAnalysis.transformations import TransformationBase +import numpy as np + +class ShiftChains(TransformationBase): + """Shift all protein chains relative to the first chain to keep them in the same box.""" + def __init__(self, prot, max_threads=1): + self.prot = prot + self.max_threads = max_threads # required by MDAnalysis + super().__init__() + + def _transform(self, ts): + chains = [seg.atoms for seg in self.prot.segments] + ref_chain = chains[0] + for chain in chains[1:]: + vec = chain.center_of_mass() - ref_chain.center_of_mass() + chain.positions -= np.rint(vec / ts.dimensions[:3]) * ts.dimensions[:3] + return ts + def make_Universe(top: pathlib.Path, trj: nc.Dataset, state: int) -> mda.Universe: """Makes a Universe and applies some transformations @@ -46,14 +65,21 @@ def make_Universe(top: pathlib.Path, trj: nc.Dataset, state: int) -> mda.Univers # - make the protein whole across periodic images between frames # - put the ligand in the closest periodic image as the protein # - align everything to minimise protein RMSD - make_whole_tr = make_whole(prot, compound="segments") + # make_whole_tr = make_whole(prot, Compound.SEGMENTS) + # Shift all chains relative to first chain to keep in same box unwrap_tr = unwrap(prot) + shift = ShiftChains(prot) + + # Make each fragment whole internally + for frag in prot.fragments: + make_whole(frag, reference_atom=frag[0]) minnie = Minimiser(prot, ligand) align = Aligner(prot) u.trajectory.add_transformations( - make_whole_tr, + # make_whole_tr, unwrap_tr, + shift, minnie, align, ) From d274b0cefcc2ad617a1590eb19f792e512e46604 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Thu, 18 Dec 2025 16:00:55 +0100 Subject: [PATCH 3/5] Add tests --- src/openfe_analysis/rmsd.py | 7 ++- src/openfe_analysis/tests/conftest.py | 8 +++ src/openfe_analysis/tests/test_rmsd.py | 82 +++++++++++++++++++++++++- 3 files changed, 92 insertions(+), 5 deletions(-) diff --git a/src/openfe_analysis/rmsd.py b/src/openfe_analysis/rmsd.py index 1d9f29e..cb52e3a 100644 --- a/src/openfe_analysis/rmsd.py +++ b/src/openfe_analysis/rmsd.py @@ -7,15 +7,13 @@ import numpy as np import tqdm from MDAnalysis.analysis import rms -from MDAnalysis.transformations import unwrap +from MDAnalysis.transformations import unwrap, TransformationBase from MDAnalysis.lib.mdamath import make_whole from numpy import typing as npt from .reader import FEReader from .transformations import Aligner, Minimiser, NoJump -from MDAnalysis.transformations import TransformationBase -import numpy as np class ShiftChains(TransformationBase): """Shift all protein chains relative to the first chain to keep them in the same box.""" @@ -149,6 +147,9 @@ def gather_rms_data( # cheeky, but we can read the PDB topology once and reuse per universe # this then only hits the PDB file once for all replicas u = make_Universe(u_top._topology, ds, state=i) + with mda.Writer("debug_after_pbc.xtc", u.atoms.n_atoms) as W: + for ts in u.trajectory[:100]: + W.write(u.atoms) prot = u.select_atoms("protein and name CA") ligand = u.select_atoms("resname UNK") diff --git a/src/openfe_analysis/tests/conftest.py b/src/openfe_analysis/tests/conftest.py index 5fa2db4..b237f47 100644 --- a/src/openfe_analysis/tests/conftest.py +++ b/src/openfe_analysis/tests/conftest.py @@ -30,6 +30,10 @@ def rbfe_skipped_data_dir() -> pathlib.Path: def simulation_nc(rbfe_output_data_dir) -> pathlib.Path: return rbfe_output_data_dir/"simulation.nc" +@pytest.fixture(scope="session") +def simulation_nc_multichain() -> pathlib.Path: + return "data/complex.nc" + @pytest.fixture(scope="session") def simulation_skipped_nc(rbfe_skipped_data_dir) -> pathlib.Path: @@ -40,6 +44,10 @@ def simulation_skipped_nc(rbfe_skipped_data_dir) -> pathlib.Path: def hybrid_system_pdb(rbfe_output_data_dir) -> pathlib.Path: return rbfe_output_data_dir/"hybrid_system.pdb" +@pytest.fixture(scope="session") +def system_pdb_multichain() -> pathlib.Path: + return "data/alchemical_system.pdb" + @pytest.fixture(scope="session") def hybrid_system_skipped_pdb(rbfe_skipped_data_dir)->pathlib.Path: diff --git a/src/openfe_analysis/tests/test_rmsd.py b/src/openfe_analysis/tests/test_rmsd.py index a108bbe..c5fcc87 100644 --- a/src/openfe_analysis/tests/test_rmsd.py +++ b/src/openfe_analysis/tests/test_rmsd.py @@ -2,8 +2,8 @@ import numpy as np import pytest from numpy.testing import assert_allclose - -from openfe_analysis.rmsd import gather_rms_data +from MDAnalysis.analysis import rms +from openfe_analysis.rmsd import gather_rms_data, make_Universe @pytest.mark.flaky(reruns=3) @@ -78,3 +78,81 @@ def test_gather_rms_data_regression_skippednc(simulation_skipped_nc, hybrid_syst [1.176307, 1.203364, 1.486987, 1.17462, 1.143457, 1.244173], rtol=1e-3, ) + +def test_multichain_com_continuity(simulation_nc_multichain, system_pdb_multichain): + u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) + prot = u.select_atoms("protein") + chains = [seg.atoms for seg in prot.segments] + assert len(chains) == 2 + + segments = prot.segments + assert len(segments) > 1, "Test requires multi-chain protein" + + chain_a = segments[0].atoms + chain_b = segments[1].atoms + + distances = [] + for ts in u.trajectory[:50]: + d = np.linalg.norm( + chain_a.center_of_mass() - chain_b.center_of_mass() + ) + distances.append(d) + + # No large frame-to-frame jumps (PBC artifacts) + jumps = np.abs(np.diff(distances)) + assert np.max(jumps) < 5.0 # Å + del u + +# def test_chain_radius_of_gyration_stable(simulation_nc_multichain, system_pdb_multichain): +# u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) +# +# protein = u.select_atoms("protein") +# chain = protein.segments[0].atoms +# +# rgs = [] +# for ts in u.trajectory[:50]: +# rgs.append(chain.radius_of_gyration()) +# +# # Chain should not explode or collapse due to PBC errors +# assert np.std(rgs) < 2.0 + +def test_rmsd_continuity(simulation_nc_multichain, system_pdb_multichain): + u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) + + prot = u.select_atoms("protein and name CA") + ref = prot.positions.copy() + + rmsds = [] + for ts in u.trajectory[:20]: + diff = prot.positions - ref + rmsd = np.sqrt((diff * diff).sum(axis=1).mean()) + rmsds.append(rmsd) + + jumps = np.abs(np.diff(rmsds)) + assert np.max(jumps) < 2.0 + del u + + +def test_rmsd_reference_is_first_frame(simulation_nc_multichain, system_pdb_multichain): + # RMS of first frame should be zero + u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) + prot = u.select_atoms("protein") + + u.trajectory[0] + ref = prot.positions.copy() + + u.trajectory[0] + r0 = rms.rmsd(prot.positions, ref, center=False, superposition=False) + + assert r0 == pytest.approx(0.0) + del u + +def test_ligand_com_continuity(simulation_nc_multichain, system_pdb_multichain): + u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) + ligand = u.select_atoms("resname UNK") + + coms = [ligand.center_of_mass() for ts in u.trajectory[:50]] + jumps = [np.linalg.norm(coms[i+1] - coms[i]) for i in range(len(coms)-1)] + + assert max(jumps) < 5.0 + del u \ No newline at end of file From f3634ddb8611c9c237823796207d20e34ddd3fdc Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Fri, 19 Dec 2025 10:25:47 +0100 Subject: [PATCH 4/5] Some fixes --- src/openfe_analysis/reader.py | 6 ++- src/openfe_analysis/rmsd.py | 11 ++--- src/openfe_analysis/tests/test_rmsd.py | 61 ++++++++++++++++---------- 3 files changed, 45 insertions(+), 33 deletions(-) diff --git a/src/openfe_analysis/reader.py b/src/openfe_analysis/reader.py index b7677bc..884550b 100644 --- a/src/openfe_analysis/reader.py +++ b/src/openfe_analysis/reader.py @@ -193,5 +193,7 @@ def _reopen(self): self._frame_index = -1 def close(self): - if self._dataset_owner: - self._dataset.close() + if self._dataset is not None: + if self._dataset_owner: + self._dataset.close() + self._dataset = None diff --git a/src/openfe_analysis/rmsd.py b/src/openfe_analysis/rmsd.py index cb52e3a..a9dbfdd 100644 --- a/src/openfe_analysis/rmsd.py +++ b/src/openfe_analysis/rmsd.py @@ -123,7 +123,6 @@ def gather_rms_data( "ligand_wander": [], "protein_2D_RMSD": [], } - ds = nc.Dataset(dataset) n_lambda = ds.dimensions["state"].size @@ -134,11 +133,11 @@ def gather_rms_data( else: n_frames = ds.dimensions["iteration"].size + if skip is None: # find skip that would give ~500 frames of output # max against 1 to avoid skip=0 case skip = max(n_frames // 500, 1) - pb = tqdm.tqdm(total=int(n_frames / skip) * n_lambda) u_top = mda.Universe(pdb_topology) @@ -147,13 +146,9 @@ def gather_rms_data( # cheeky, but we can read the PDB topology once and reuse per universe # this then only hits the PDB file once for all replicas u = make_Universe(u_top._topology, ds, state=i) - with mda.Writer("debug_after_pbc.xtc", u.atoms.n_atoms) as W: - for ts in u.trajectory[:100]: - W.write(u.atoms) prot = u.select_atoms("protein and name CA") ligand = u.select_atoms("resname UNK") - # save coordinates for 2D RMSD matrix # TODO: Some smart guard to avoid allocating a silly amount of memory? prot2d = np.empty((len(u.trajectory[::skip]), len(prot), 3), dtype=np.float32) @@ -172,6 +167,7 @@ def gather_rms_data( pb.update() if prot: + # prot2d[ts_i] = prot.positions prot2d[ts_i, :, :] = prot.positions this_protein_rmsd.append( rms.rmsd( @@ -197,6 +193,7 @@ def gather_rms_data( # ignores PBC, but we've already centered the traj mda.lib.distances.calc_bonds(ligand.center_of_mass(), ligand_initial_com) ) + # ts_i += 1 if prot: # can ignore weights here as it's all Ca @@ -208,7 +205,7 @@ def gather_rms_data( output["ligand_wander"].append(this_ligand_wander) output["time(ps)"] = list(np.arange(len(u.trajectory))[::skip] * u.trajectory.dt) - + ds.close() return output diff --git a/src/openfe_analysis/tests/test_rmsd.py b/src/openfe_analysis/tests/test_rmsd.py index c5fcc87..c284d85 100644 --- a/src/openfe_analysis/tests/test_rmsd.py +++ b/src/openfe_analysis/tests/test_rmsd.py @@ -1,12 +1,29 @@ import netCDF4 as nc import numpy as np import pytest +from itertools import islice from numpy.testing import assert_allclose from MDAnalysis.analysis import rms from openfe_analysis.rmsd import gather_rms_data, make_Universe +@pytest.fixture +def mda_universe(system_pdb_multichain, simulation_nc_multichain): + """ + Safely create and destroy an MDAnalysis Universe. + + Guarantees: + - NetCDF file is opened exactly once + """ + u = make_Universe( + system_pdb_multichain, + simulation_nc_multichain, + state=0, + ) + + yield u + -@pytest.mark.flaky(reruns=3) +@pytest.mark.flaky(reruns=1) def test_gather_rms_data_regression(simulation_nc, hybrid_system_pdb): output = gather_rms_data( hybrid_system_pdb, @@ -43,7 +60,7 @@ def test_gather_rms_data_regression(simulation_nc, hybrid_system_pdb): ) -@pytest.mark.flaky(reruns=3) +@pytest.mark.flaky(reruns=1) def test_gather_rms_data_regression_skippednc(simulation_skipped_nc, hybrid_system_skipped_pdb): output = gather_rms_data( hybrid_system_skipped_pdb, @@ -79,8 +96,8 @@ def test_gather_rms_data_regression_skippednc(simulation_skipped_nc, hybrid_syst rtol=1e-3, ) -def test_multichain_com_continuity(simulation_nc_multichain, system_pdb_multichain): - u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) +def test_multichain_com_continuity(mda_universe): + u = mda_universe prot = u.select_atoms("protein") chains = [seg.atoms for seg in prot.segments] assert len(chains) == 2 @@ -92,7 +109,7 @@ def test_multichain_com_continuity(simulation_nc_multichain, system_pdb_multicha chain_b = segments[1].atoms distances = [] - for ts in u.trajectory[:50]: + for ts in islice(u.trajectory, 20): d = np.linalg.norm( chain_a.center_of_mass() - chain_b.center_of_mass() ) @@ -101,7 +118,7 @@ def test_multichain_com_continuity(simulation_nc_multichain, system_pdb_multicha # No large frame-to-frame jumps (PBC artifacts) jumps = np.abs(np.diff(distances)) assert np.max(jumps) < 5.0 # Å - del u + u.trajectory.close() # def test_chain_radius_of_gyration_stable(simulation_nc_multichain, system_pdb_multichain): # u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) @@ -116,43 +133,39 @@ def test_multichain_com_continuity(simulation_nc_multichain, system_pdb_multicha # # Chain should not explode or collapse due to PBC errors # assert np.std(rgs) < 2.0 -def test_rmsd_continuity(simulation_nc_multichain, system_pdb_multichain): - u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) +def test_rmsd_continuity(mda_universe): + u = mda_universe prot = u.select_atoms("protein and name CA") ref = prot.positions.copy() rmsds = [] - for ts in u.trajectory[:20]: + for ts in islice(u.trajectory, 20): diff = prot.positions - ref rmsd = np.sqrt((diff * diff).sum(axis=1).mean()) rmsds.append(rmsd) jumps = np.abs(np.diff(rmsds)) assert np.max(jumps) < 2.0 - del u + u.trajectory.close() - -def test_rmsd_reference_is_first_frame(simulation_nc_multichain, system_pdb_multichain): - # RMS of first frame should be zero - u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) +def test_rmsd_reference_is_first_frame(mda_universe): + u = mda_universe prot = u.select_atoms("protein") - u.trajectory[0] + ts = next(iter(u.trajectory)) # SAFE ref = prot.positions.copy() - u.trajectory[0] - r0 = rms.rmsd(prot.positions, ref, center=False, superposition=False) - - assert r0 == pytest.approx(0.0) - del u + rmsd = np.sqrt(((prot.positions - ref) ** 2).mean()) + assert rmsd == 0.0 + u.trajectory.close() -def test_ligand_com_continuity(simulation_nc_multichain, system_pdb_multichain): - u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) +def test_ligand_com_continuity(mda_universe): + u = mda_universe ligand = u.select_atoms("resname UNK") - coms = [ligand.center_of_mass() for ts in u.trajectory[:50]] + coms = [ligand.center_of_mass() for ts in islice(u.trajectory, 20)] jumps = [np.linalg.norm(coms[i+1] - coms[i]) for i in range(len(coms)-1)] assert max(jumps) < 5.0 - del u \ No newline at end of file + u.trajectory.close() From e92adb3286769458db77930ed19b37f0879a3653 Mon Sep 17 00:00:00 2001 From: hannahbaumann Date: Fri, 19 Dec 2025 11:04:16 +0100 Subject: [PATCH 5/5] Add another test --- src/openfe_analysis/rmsd.py | 8 +++----- src/openfe_analysis/tests/test_rmsd.py | 25 +++++++++++++------------ 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/src/openfe_analysis/rmsd.py b/src/openfe_analysis/rmsd.py index a9dbfdd..cb4cff3 100644 --- a/src/openfe_analysis/rmsd.py +++ b/src/openfe_analysis/rmsd.py @@ -63,7 +63,6 @@ def make_Universe(top: pathlib.Path, trj: nc.Dataset, state: int) -> mda.Univers # - make the protein whole across periodic images between frames # - put the ligand in the closest periodic image as the protein # - align everything to minimise protein RMSD - # make_whole_tr = make_whole(prot, Compound.SEGMENTS) # Shift all chains relative to first chain to keep in same box unwrap_tr = unwrap(prot) shift = ShiftChains(prot) @@ -75,7 +74,6 @@ def make_Universe(top: pathlib.Path, trj: nc.Dataset, state: int) -> mda.Univers align = Aligner(prot) u.trajectory.add_transformations( - # make_whole_tr, unwrap_tr, shift, minnie, @@ -123,6 +121,7 @@ def gather_rms_data( "ligand_wander": [], "protein_2D_RMSD": [], } + ds = nc.Dataset(dataset) n_lambda = ds.dimensions["state"].size @@ -133,11 +132,11 @@ def gather_rms_data( else: n_frames = ds.dimensions["iteration"].size - if skip is None: # find skip that would give ~500 frames of output # max against 1 to avoid skip=0 case skip = max(n_frames // 500, 1) + pb = tqdm.tqdm(total=int(n_frames / skip) * n_lambda) u_top = mda.Universe(pdb_topology) @@ -149,6 +148,7 @@ def gather_rms_data( prot = u.select_atoms("protein and name CA") ligand = u.select_atoms("resname UNK") + # save coordinates for 2D RMSD matrix # TODO: Some smart guard to avoid allocating a silly amount of memory? prot2d = np.empty((len(u.trajectory[::skip]), len(prot), 3), dtype=np.float32) @@ -167,7 +167,6 @@ def gather_rms_data( pb.update() if prot: - # prot2d[ts_i] = prot.positions prot2d[ts_i, :, :] = prot.positions this_protein_rmsd.append( rms.rmsd( @@ -193,7 +192,6 @@ def gather_rms_data( # ignores PBC, but we've already centered the traj mda.lib.distances.calc_bonds(ligand.center_of_mass(), ligand_initial_com) ) - # ts_i += 1 if prot: # can ignore weights here as it's all Ca diff --git a/src/openfe_analysis/tests/test_rmsd.py b/src/openfe_analysis/tests/test_rmsd.py index c284d85..898d44d 100644 --- a/src/openfe_analysis/tests/test_rmsd.py +++ b/src/openfe_analysis/tests/test_rmsd.py @@ -120,18 +120,19 @@ def test_multichain_com_continuity(mda_universe): assert np.max(jumps) < 5.0 # Å u.trajectory.close() -# def test_chain_radius_of_gyration_stable(simulation_nc_multichain, system_pdb_multichain): -# u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) -# -# protein = u.select_atoms("protein") -# chain = protein.segments[0].atoms -# -# rgs = [] -# for ts in u.trajectory[:50]: -# rgs.append(chain.radius_of_gyration()) -# -# # Chain should not explode or collapse due to PBC errors -# assert np.std(rgs) < 2.0 +def test_chain_radius_of_gyration_stable(simulation_nc_multichain, system_pdb_multichain): + u = make_Universe(system_pdb_multichain, simulation_nc_multichain, state=0) + + protein = u.select_atoms("protein") + chain = protein.segments[0].atoms + + rgs = [] + for ts in u.trajectory[:50]: + rgs.append(chain.radius_of_gyration()) + + # Chain should not explode or collapse due to PBC errors + assert np.std(rgs) < 2.0 + u.trajectory.close() def test_rmsd_continuity(mda_universe): u = mda_universe