From a73586fdaa94b12113c27e835e8872cf1cd7962f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Mon, 10 Jul 2023 17:52:32 -0400 Subject: [PATCH 1/8] Initial and final topologies serialized per phase. --- perses/app/relative_setup.py | 74 +++++++++++++++++--- perses/tests/test_relative_setup.py | 100 +++++++++++++++------------- 2 files changed, 120 insertions(+), 54 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index a830e9986..c54471ca8 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -587,6 +587,9 @@ def __init__(self, self._vacuum_reverse_neglected_angles = self._geometry_engine.reverse_neglected_angle_terms self._vacuum_geometry_engine = copy.deepcopy(self._geometry_engine) + # Serialize topologies and positions in PDB files + self._serialize_topologies() + def _setup_complex_phase(self): """ Runs setup on the protein/receptor file for relative simulations @@ -610,6 +613,8 @@ def _setup_complex_phase(self): raise ValueError("You need to provide either a protein pdb or a receptor mol2 to run a complex simulation.") self._complex_md_topology_old = self._receptor_md_topology_old.join(self._ligand_md_topology_old) + # receptor is the same for "old" and "new" -- only for small molecule pipeline + self._complex_md_topology_new = self._receptor_md_topology_old.join(self._ligand_md_topology_new) n_atoms_spectators = 0 if self._spectator_filenames: @@ -618,14 +623,22 @@ def _setup_complex_phase(self): self._complex_md_topology_old = self._complex_md_topology_old.join(spectator_topology) n_atoms_spectators += spectator_topology.n_atoms self._complex_topology_old = self._complex_md_topology_old.to_openmm() + self._complex_topology_new = self._complex_md_topology_new.to_openmm() n_atoms_total_old = self._complex_topology_old.getNumAtoms() + n_atoms_total_new = self._complex_topology_new.getNumAtoms() n_atoms_protein_old = self._receptor_topology_old.getNumAtoms() + n_atoms_protein_new = n_atoms_protein_old # Assumes it's the same target n_atoms_ligand_old = n_atoms_total_old - n_atoms_protein_old - n_atoms_spectators + n_atoms_ligand_new = n_atoms_total_new - n_atoms_protein_new - n_atoms_spectators self._complex_positions_old = unit.Quantity(np.zeros([n_atoms_total_old, 3]), unit=unit.nanometers) self._complex_positions_old[:n_atoms_protein_old, :] = self._receptor_positions_old self._complex_positions_old[n_atoms_protein_old:n_atoms_protein_old+n_atoms_ligand_old, :] = self._ligand_positions_old + self._complex_positions_new = unit.Quantity(np.zeros([n_atoms_total_new, 3]), unit=unit.nanometers) + self._complex_positions_new[:n_atoms_protein_new, :] = self._receptor_positions_old # Assumes it's the same target + self._complex_positions_new[n_atoms_protein_new:n_atoms_protein_new + n_atoms_ligand_new, + :] = self._ligand_positions_new if self._spectator_filenames: start = n_atoms_protein_old+n_atoms_ligand_old @@ -674,6 +687,9 @@ def _generate_solvent_topologies(self, topology_proposal, old_positions): old_solvated_topology, old_solvated_positions, old_solvated_system = self._solvate_system( old_ligand_topology.to_openmm(), old_ligand_positions, phase='solvent', box_dimensions=self._solvent_box_dimensions, model=self._solvent_model) + # Update attributes + self._ligand_topology_old_solvated = old_solvated_topology + self._ligand_positions_old_solvated = old_solvated_positions old_solvated_md_topology = md.Topology.from_openmm(old_solvated_topology) @@ -847,7 +863,6 @@ def _solvate_system(self, topology, positions, model='tip3p',phase='complex', bo The parameterized system, containing a barostat if one was specified. """ # DEBUG: Write PDB file being fed into Modeller to check why MOL isn't being matched - from simtk.openmm.app import PDBFile modeller = app.Modeller(topology, positions) # retaining protein protonation from input files #hs = [atom for atom in modeller.topology.atoms() if atom.element.symbol in ['H'] and atom.residue.name not in ['MOL','OLD','NEW']] @@ -881,13 +896,6 @@ def _solvate_system(self, topology, positions, model='tip3p',phase='complex', bo solvated_system = self._system_generator.create_system(solvated_topology) _logger.info(f"\tSystem parameterized") - if self._trajectory_directory is not None and self._trajectory_prefix is not None: - pdb_filename = AnyPath(f"{self._trajectory_directory}/{self._trajectory_prefix}-{phase}.pdb") - with open(pdb_filename, 'w') as outfile: - PDBFile.writeFile(solvated_topology, solvated_positions, outfile) - else: - _logger.info('Both trajectory_directory and trajectory_prefix need to be provided to save .pdb') - return solvated_topology, solvated_positions, solvated_system def _handle_charge_changes(self, topology_proposal, new_positions): @@ -922,6 +930,56 @@ def _handle_charge_changes(self, topology_proposal, new_positions): # modify the topology proposal modify_atom_classes(new_water_indices_to_ionize, topology_proposal) + def _serialize_topologies(self): + """ + Serializes topologies and positions as PDB files. + + Notes + ----- + - The topologies and positions are serialized in the specified directory using the specified prefix. + - The directory will be created if it doesn't already exist. + - The topologies and positions are stored as PDB files in the "models" subdirectory of the trajectory directory. + - The following topologies and positions are serialized: + - "solvent_old": old solvent topology and positions + - "solvent_old_solvated": old solvated solvent topology and positions + - "solvent_new": new solvent topology and positions + - "complex_old": old complex topology and positions + - "complex_old_solvated": old solvated complex topology and positions + - "complex_new": new complex topology and positions + - Both the trajectory directory and trajectory prefix need to be provided in order to serialize the topologies. + - If any of the directories or files already exist, they will be overwritten. + + Usage + ----- + To serialize the topologies and positions, make sure to set the `trajectory_directory` and `trajectory_prefix` + attributes before calling this method. + """ + from openmm.app import PDBFile + if self._trajectory_directory is not None and self._trajectory_prefix is not None: + topologies_to_serialize = { + "solvent_old": {"topology": self._ligand_topology_old, "positions": self._ligand_positions_old}, + "solvent_old_solvated": {"topology": self._ligand_topology_old_solvated, + "positions": self._ligand_positions_old_solvated}, + "solvent_new": {"topology": self._ligand_topology_new, "positions": self._ligand_positions_new}, + "complex_old": {"topology": self._complex_topology_old, "positions": self._complex_positions_old}, + "complex_old_solvated": {"topology": self._complex_topology_old_solvated, + "positions": self._complex_positions_old_solvated}, + "complex_new": {"topology": self._complex_topology_new, "positions": self._complex_positions_new}} + # Serialize in "models" subdir -- create if doesn't exist + models_path = f"{self._trajectory_directory}/models" + if not os.path.exists(models_path): + os.makedirs(models_path) + for phase, phase_data in topologies_to_serialize.items(): + pdb_filename = AnyPath(f"{models_path}/{self._trajectory_prefix}_{phase}.pdb") + with open(pdb_filename, 'w') as outfile: + PDBFile.writeFile(phase_data["topology"], phase_data["positions"], outfile) + else: + _logger.info( + 'Not serializing topologies. Both trajectory_directory and trajectory_prefix need to be provided to ' + 'store topology .pdb files') + + + @property def complex_topology_proposal(self): return self._complex_topology_proposal diff --git a/perses/tests/test_relative_setup.py b/perses/tests/test_relative_setup.py index b1029bedb..5bdcabcaa 100644 --- a/perses/tests/test_relative_setup.py +++ b/perses/tests/test_relative_setup.py @@ -7,7 +7,6 @@ from perses.app import setup_relative_calculation import mdtraj as md from openmmtools import states, alchemy, testsystems, cache -from unittest import skipIf import pytest from perses.tests.utils import enter_temp_directory @@ -19,7 +18,6 @@ 'lambda_electrostatics' : 'lambda', } - def generate_example_waterbox_states(temperature=300.0*unit.kelvin, pressure=1.0*unit.atmosphere): """ This is a convenience function to generate a CompoundThermodynamicState and SamplerState to use in other tests. @@ -260,47 +258,6 @@ def warn(*args, **kwargs): found_tla = True assert found_tla == True, 'Spectator TLA not in old topology' - -def test_host_guest_deterministic_geometries(): - """ - execute the `RelativeFEPSetup` with geometries specified for a host-guest relative free energy pair - """ - from perses.app.relative_setup import RelativeFEPSetup - - # Setup directory - ligand_sdf = resource_filename("perses", "data/given-geometries/ligands.sdf") - host_pdb = resource_filename("perses", "data/given-geometries/receptor.pdb") - - setup = RelativeFEPSetup( - ligand_input = ligand_sdf, - old_ligand_index=0, - new_ligand_index=1, - forcefield_files = ['amber/ff14SB.xml','amber/tip3p_standard.xml','amber/tip3p_HFE_multivalent.xml'], - phases = ['complex', 'solvent', 'vacuum'], - protein_pdb_filename=host_pdb, - receptor_mol2_filename=None, - pressure=1.0 * unit.atmosphere, - temperature=300.0 * unit.kelvin, - solvent_padding=9.0 * unit.angstroms, - ionic_strength=0.15 * unit.molar, - hmass=3*unit.amus, - neglect_angles=False, - map_strength='default', - atom_expr=None, - bond_expr=None, - anneal_14s=False, - small_molecule_forcefield='gaff-2.11', - small_molecule_parameters_cache=None, - trajectory_directory=None, - trajectory_prefix=None, - spectator_filenames=None, - nonbonded_method = 'PME', - complex_box_dimensions=None, - solvent_box_dimensions=None, - remove_constraints=False, - use_given_geometries = True - ) - def test_relative_setup_charge_change(): """ execute `RelativeFEPSetup` in solvent/complex phase on a charge change and assert that the modified new system and old system charge difference is zero. @@ -357,7 +314,6 @@ def test_relative_setup_solvent_padding(): Check that the user inputted solvent_padding argument to `RelativeFEPSetup` actually changes the padding for the solvent phase """ from perses.app.relative_setup import RelativeFEPSetup - import numpy as np input_solvent_padding = 1.7 * unit.nanometers smiles_filename = resource_filename("perses", os.path.join("data", "test.smi")) @@ -409,6 +365,58 @@ def test_relative_setup_list_ligand_input(): assert ligand_input_size == expected_input_size, f"There should be {expected_input_size} ligand input files, " \ f"receiving {ligand_input_size} files." +def test_relative_setup_topologies_serialization(tmp_path): + """ + Check that topologies + """ + from perses.app.relative_setup import RelativeFEPSetup -# if __name__=="__main__": -# test_run_cdk2_iterations_repex() + # Setup directory + ligand_sdf = resource_filename("perses", "data/given-geometries/ligands.sdf") + host_pdb = resource_filename("perses", "data/given-geometries/receptor.pdb") + + fe_setup = RelativeFEPSetup( + ligand_input=ligand_sdf, + old_ligand_index=0, + new_ligand_index=1, + forcefield_files=['amber/ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml'], + phases=['complex', 'solvent', 'vacuum'], + protein_pdb_filename=host_pdb, + receptor_mol2_filename=None, + pressure=1.0 * unit.atmosphere, + temperature=300.0 * unit.kelvin, + solvent_padding=9.0 * unit.angstroms, + ionic_strength=0.15 * unit.molar, + hmass=3 * unit.amus, + neglect_angles=False, + map_strength='default', + atom_expr=None, + bond_expr=None, + anneal_14s=False, + small_molecule_forcefield='openff-2.1.0', + small_molecule_parameters_cache=None, + trajectory_directory=tmp_path, + trajectory_prefix="host-guest", + spectator_filenames=None, + nonbonded_method='PME', + complex_box_dimensions=None, + solvent_box_dimensions=None, + remove_constraints=False, + use_given_geometries=True + ) + + expected_files = [ + 'host-guest_solvent_old.pdb', + 'host-guest_solvent_old_solvated.pdb', + 'host-guest_solvent_new.pdb', + 'host-guest_complex_old.pdb', + 'host-guest_complex_old_solvated.pdb', + 'host-guest_complex_new.pdb' + ] + for file in expected_files: + pdb_filename = f"{tmp_path}/models/{file}" + assert os.path.exists(pdb_filename) + + # Assert that no extra files are created + files_in_directory = os.listdir(f"{tmp_path}/models") + assert len(files_in_directory) == len(expected_files) From 9ada08e5a8ffe8a767bf4fc8290045f4ecd1bc0d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Thu, 27 Jul 2023 14:47:18 -0400 Subject: [PATCH 2/8] Using properties instead of private attrs. --- perses/app/relative_setup.py | 40 +++++++++++++++++------------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index c54471ca8..5d5e4b80f 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -587,8 +587,8 @@ def __init__(self, self._vacuum_reverse_neglected_angles = self._geometry_engine.reverse_neglected_angle_terms self._vacuum_geometry_engine = copy.deepcopy(self._geometry_engine) - # Serialize topologies and positions in PDB files - self._serialize_topologies() + # Store topologies and positions in PDB files + self._store_topologies() def _setup_complex_phase(self): """ @@ -930,41 +930,39 @@ def _handle_charge_changes(self, topology_proposal, new_positions): # modify the topology proposal modify_atom_classes(new_water_indices_to_ionize, topology_proposal) - def _serialize_topologies(self): + def _store_topologies(self): """ - Serializes topologies and positions as PDB files. + Stores solvated topologies and positions in PDB files. Notes ----- - - The topologies and positions are serialized in the specified directory using the specified prefix. + - The topologies and positions are stored in the specified directory using the global trajectory prefix. - The directory will be created if it doesn't already exist. - The topologies and positions are stored as PDB files in the "models" subdirectory of the trajectory directory. - - The following topologies and positions are serialized: + - The following solvated topologies and positions are stored: - "solvent_old": old solvent topology and positions - - "solvent_old_solvated": old solvated solvent topology and positions - "solvent_new": new solvent topology and positions - "complex_old": old complex topology and positions - - "complex_old_solvated": old solvated complex topology and positions - "complex_new": new complex topology and positions - - Both the trajectory directory and trajectory prefix need to be provided in order to serialize the topologies. + - Both the trajectory directory and trajectory prefix need to be provided in order to store the topologies. - If any of the directories or files already exist, they will be overwritten. Usage ----- - To serialize the topologies and positions, make sure to set the `trajectory_directory` and `trajectory_prefix` + To store the topologies and positions, make sure to set the `trajectory_directory` and `trajectory_prefix` attributes before calling this method. """ from openmm.app import PDBFile if self._trajectory_directory is not None and self._trajectory_prefix is not None: topologies_to_serialize = { - "solvent_old": {"topology": self._ligand_topology_old, "positions": self._ligand_positions_old}, - "solvent_old_solvated": {"topology": self._ligand_topology_old_solvated, - "positions": self._ligand_positions_old_solvated}, - "solvent_new": {"topology": self._ligand_topology_new, "positions": self._ligand_positions_new}, - "complex_old": {"topology": self._complex_topology_old, "positions": self._complex_positions_old}, - "complex_old_solvated": {"topology": self._complex_topology_old_solvated, - "positions": self._complex_positions_old_solvated}, - "complex_new": {"topology": self._complex_topology_new, "positions": self._complex_positions_new}} + "solvent_old": {"topology": self.solvent_topology_proposal.old_topology, + "positions": self.solvent_old_positions}, + "solvent_new": {"topology": self.solvent_topology_proposal.new_topology, + "positions": self.solvent_new_positions}, + "complex_old": {"topology": self.complex_topology_proposal.old_topology, + "positions": self.complex_old_positions}, + "complex_new": {"topology": self.complex_topology_proposal.new_topology, + "positions": self.complex_new_positions}} # Serialize in "models" subdir -- create if doesn't exist models_path = f"{self._trajectory_directory}/models" if not os.path.exists(models_path): @@ -974,9 +972,9 @@ def _serialize_topologies(self): with open(pdb_filename, 'w') as outfile: PDBFile.writeFile(phase_data["topology"], phase_data["positions"], outfile) else: - _logger.info( - 'Not serializing topologies. Both trajectory_directory and trajectory_prefix need to be provided to ' - 'store topology .pdb files') + _logger.warning( + 'Not storing topologies. Both trajectory_directory and trajectory_prefix arguments need to' + ' be provided to store topology .pdb files.') From 36a74582182471d879f66935f61ec610457f1aed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Thu, 27 Jul 2023 15:26:59 -0400 Subject: [PATCH 3/8] Fix test --- perses/tests/test_relative_setup.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/perses/tests/test_relative_setup.py b/perses/tests/test_relative_setup.py index 5bdcabcaa..7a3e8e0e8 100644 --- a/perses/tests/test_relative_setup.py +++ b/perses/tests/test_relative_setup.py @@ -365,9 +365,9 @@ def test_relative_setup_list_ligand_input(): assert ligand_input_size == expected_input_size, f"There should be {expected_input_size} ligand input files, " \ f"receiving {ligand_input_size} files." -def test_relative_setup_topologies_serialization(tmp_path): +def test_relative_setup_topologies_storage(tmp_path): """ - Check that topologies + Check that topologies get stored """ from perses.app.relative_setup import RelativeFEPSetup @@ -407,10 +407,8 @@ def test_relative_setup_topologies_serialization(tmp_path): expected_files = [ 'host-guest_solvent_old.pdb', - 'host-guest_solvent_old_solvated.pdb', 'host-guest_solvent_new.pdb', 'host-guest_complex_old.pdb', - 'host-guest_complex_old_solvated.pdb', 'host-guest_complex_new.pdb' ] for file in expected_files: From 2e7aa04d8445527e11f050147f0e7fb529d15ed7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Thu, 27 Jul 2023 15:35:14 -0400 Subject: [PATCH 4/8] Remove uneeded code/attributes --- perses/app/relative_setup.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index 021965ed5..ab153b764 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -613,8 +613,6 @@ def _setup_complex_phase(self): raise ValueError("You need to provide either a protein pdb or a receptor mol2 to run a complex simulation.") self._complex_md_topology_old = self._receptor_md_topology_old.join(self._ligand_md_topology_old) - # receptor is the same for "old" and "new" -- only for small molecule pipeline - self._complex_md_topology_new = self._receptor_md_topology_old.join(self._ligand_md_topology_new) n_atoms_spectators = 0 if self._spectator_filenames: @@ -623,22 +621,14 @@ def _setup_complex_phase(self): self._complex_md_topology_old = self._complex_md_topology_old.join(spectator_topology) n_atoms_spectators += spectator_topology.n_atoms self._complex_topology_old = self._complex_md_topology_old.to_openmm() - self._complex_topology_new = self._complex_md_topology_new.to_openmm() n_atoms_total_old = self._complex_topology_old.getNumAtoms() - n_atoms_total_new = self._complex_topology_new.getNumAtoms() n_atoms_protein_old = self._receptor_topology_old.getNumAtoms() - n_atoms_protein_new = n_atoms_protein_old # Assumes it's the same target n_atoms_ligand_old = n_atoms_total_old - n_atoms_protein_old - n_atoms_spectators - n_atoms_ligand_new = n_atoms_total_new - n_atoms_protein_new - n_atoms_spectators self._complex_positions_old = unit.Quantity(np.zeros([n_atoms_total_old, 3]), unit=unit.nanometers) self._complex_positions_old[:n_atoms_protein_old, :] = self._receptor_positions_old self._complex_positions_old[n_atoms_protein_old:n_atoms_protein_old+n_atoms_ligand_old, :] = self._ligand_positions_old - self._complex_positions_new = unit.Quantity(np.zeros([n_atoms_total_new, 3]), unit=unit.nanometers) - self._complex_positions_new[:n_atoms_protein_new, :] = self._receptor_positions_old # Assumes it's the same target - self._complex_positions_new[n_atoms_protein_new:n_atoms_protein_new + n_atoms_ligand_new, - :] = self._ligand_positions_new if self._spectator_filenames: start = n_atoms_protein_old+n_atoms_ligand_old @@ -687,9 +677,6 @@ def _generate_solvent_topologies(self, topology_proposal, old_positions): old_solvated_topology, old_solvated_positions, old_solvated_system = self._solvate_system( old_ligand_topology.to_openmm(), old_ligand_positions, phase='solvent', box_dimensions=self._solvent_box_dimensions, model=self._solvent_model) - # Update attributes - self._ligand_topology_old_solvated = old_solvated_topology - self._ligand_positions_old_solvated = old_solvated_positions old_solvated_md_topology = md.Topology.from_openmm(old_solvated_topology) From be7f52a74d61000fee7e9959e759ca54c0fcfb5d Mon Sep 17 00:00:00 2001 From: Mike Henry <11765982+mikemhenry@users.noreply.github.com> Date: Thu, 27 Jul 2023 13:03:51 -0700 Subject: [PATCH 5/8] bump ci From f0b23bd6614ad19513d07b189239823ceec99d6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Thu, 27 Jul 2023 17:34:06 -0400 Subject: [PATCH 6/8] Store phase topologies separately --- perses/app/relative_setup.py | 57 +++++++++++++++++++++--------------- 1 file changed, 34 insertions(+), 23 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index ab153b764..012a3933e 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -461,6 +461,9 @@ def __init__(self, else: _logger.info("Skipping counterion") + # Store phase topology and positions in PDB file + self._store_phase_topologies(phase="complex") + if 'solvent' in phases: _logger.info(f"Detected solvent...") @@ -521,6 +524,9 @@ def __init__(self, else: _logger.info("Skipping counterion") + # Store phase topology and positions in PDB file + self._store_phase_topologies(phase="solvent") + if 'vacuum' in phases: _logger.info(f"Detected solvent...") # need to change nonbonded cutoff and remove barostat for vacuum leg @@ -586,9 +592,8 @@ def __init__(self, self._vacuum_forward_neglected_angles = self._geometry_engine.forward_neglected_angle_terms self._vacuum_reverse_neglected_angles = self._geometry_engine.reverse_neglected_angle_terms self._vacuum_geometry_engine = copy.deepcopy(self._geometry_engine) - - # Store topologies and positions in PDB files - self._store_topologies() + # Store phase topology and positions in PDB file + self._store_phase_topologies(phase="vacuum") def _setup_complex_phase(self): """ @@ -917,22 +922,18 @@ def _handle_charge_changes(self, topology_proposal, new_positions): # modify the topology proposal modify_atom_classes(new_water_indices_to_ionize, topology_proposal) - def _store_topologies(self): + def _store_phase_topologies(self, phase: str): """ - Stores solvated topologies and positions in PDB files. + Stores topologies and positions in PDB file for the given the phase. Stores both solvent and solute + whenever possible (only "solute" in vacuum). Notes ----- - The topologies and positions are stored in the specified directory using the global trajectory prefix. - The directory will be created if it doesn't already exist. - The topologies and positions are stored as PDB files in the "models" subdirectory of the trajectory directory. - - The following solvated topologies and positions are stored: - - "solvent_old": old solvent topology and positions - - "solvent_new": new solvent topology and positions - - "complex_old": old complex topology and positions - - "complex_new": new complex topology and positions - Both the trajectory directory and trajectory prefix need to be provided in order to store the topologies. - - If any of the directories or files already exist, they will be overwritten. + - If any of the directories or files already exist, they will get overwritten. Usage ----- @@ -941,21 +942,31 @@ def _store_topologies(self): """ from openmm.app import PDBFile if self._trajectory_directory is not None and self._trajectory_prefix is not None: - topologies_to_serialize = { - "solvent_old": {"topology": self.solvent_topology_proposal.old_topology, - "positions": self.solvent_old_positions}, - "solvent_new": {"topology": self.solvent_topology_proposal.new_topology, - "positions": self.solvent_new_positions}, - "complex_old": {"topology": self.complex_topology_proposal.old_topology, - "positions": self.complex_old_positions}, - "complex_new": {"topology": self.complex_topology_proposal.new_topology, - "positions": self.complex_new_positions}} - # Serialize in "models" subdir -- create if doesn't exist + if phase == "complex": + topology_proposal = self.complex_topology_proposal + old_positions = self.complex_old_positions + new_positions = self.complex_new_positions + elif phase == "solvent": + topology_proposal = self.solvent_topology_proposal + old_positions = self.solvent_old_positions + new_positions = self.solvent_new_positions + elif phase == "vacuum": + topology_proposal = self.vacuum_topology_proposal + old_positions = self.vacuum_old_positions + new_positions = self.vacuum_new_positions + + topologies_to_store = { + f"{phase}_old": {"topology": topology_proposal.old_topology, + "positions": old_positions}, + f"{phase}_new": {"topology": topology_proposal.new_topology, + "positions": new_positions}, + } + # Store in "models" subdir -- create if doesn't exist models_path = f"{self._trajectory_directory}/models" if not os.path.exists(models_path): os.makedirs(models_path) - for phase, phase_data in topologies_to_serialize.items(): - pdb_filename = AnyPath(f"{models_path}/{self._trajectory_prefix}_{phase}.pdb") + for phase_key, phase_data in topologies_to_store.items(): + pdb_filename = AnyPath(f"{models_path}/{self._trajectory_prefix}_{phase_key}.pdb") with open(pdb_filename, 'w') as outfile: PDBFile.writeFile(phase_data["topology"], phase_data["positions"], outfile) else: From 4a593bee87433202455b1065c98029a5f1364217 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Thu, 27 Jul 2023 17:34:20 -0400 Subject: [PATCH 7/8] Fix tests. vacuum topologies expected. --- perses/tests/test_relative_setup.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/perses/tests/test_relative_setup.py b/perses/tests/test_relative_setup.py index 7a3e8e0e8..da53fac92 100644 --- a/perses/tests/test_relative_setup.py +++ b/perses/tests/test_relative_setup.py @@ -406,10 +406,12 @@ def test_relative_setup_topologies_storage(tmp_path): ) expected_files = [ + 'host-guest_vacuum_old.pdb', + 'host-guest_vacuum_new.pdb', 'host-guest_solvent_old.pdb', 'host-guest_solvent_new.pdb', 'host-guest_complex_old.pdb', - 'host-guest_complex_new.pdb' + 'host-guest_complex_new.pdb', ] for file in expected_files: pdb_filename = f"{tmp_path}/models/{file}" From 499b39ea3b61d4b21377d65343bd7afbe10517bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Iv=C3=A1n=20Pulido?= <2949729+ijpulidos@users.noreply.github.com> Date: Thu, 27 Jul 2023 17:48:16 -0400 Subject: [PATCH 8/8] Better docstring. --- perses/app/relative_setup.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index 012a3933e..e5bba621a 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -924,8 +924,10 @@ def _handle_charge_changes(self, topology_proposal, new_positions): def _store_phase_topologies(self, phase: str): """ - Stores topologies and positions in PDB file for the given the phase. Stores both solvent and solute - whenever possible (only "solute" in vacuum). + Stores topologies and positions in PDB file for the given the phase, for both the initial (old) and final (new) + states. Stores both solvent and solute whenever possible (only "solute" in vacuum). + + Generates two PDB files, one for each state (old/new). Notes -----